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PREFACE 

The overall objective of this study was to test processing tech- 
niques for deducing cloud motion vectors from overlapped portions of 
pairs of pictures made from meteorological satellites. 

The approach was to examine some previously developed comDuter 
programs, test them with artificial image pairs, evaluate their methods 
and prepare a procedure for deducing pattern motion. Experimental data 
would be selected from polar orbiting satellite exoeriments and routines 
developed to prepare that data in format for input to the motion * 
computation procedure. Actual data samples would then be used to 
evaluate the procedure and it would be modified accordingly. 

Highlights incorporated into the program were: 

a. Suppression of spurious cross covariance products by bounding 
the window area with neutral values. 

b. Replacement of missing data with noise. 

c. Normalization of cross correlations according to scale. 

d. Exclusion of data which while valid to the sensor constitute 
noise to the cloud motion computation. 

The completed program was applied to the selected data samples with 
different parametric values to measure its derivation of motion estimates, 
timing requirements, and its limitations. Vector computations were 
compared to ground truth data obtained independently as well as to 
manually estimated motions from pseudo color photographic images. 
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Section 1 
PROJECT HISTORY 

The Investigation of means of deriving quantitative measures of 
atmospheric motion by IBM scientists began in 1964-1965 with ESSA Contract 
CWb 1098, entitled "Wind Shear Measurement by Satellite Cloud Tracking." 
This study, performed by N. Woodrick, L. Bodin, and J. Leese, analyzed the 
problem of obtaining vertical wind shear information from satellite 
observations of the motions of clouds relative to each other. 

Subsequently, during 1967 and early 1968, two studies were internally 
funded by the IBM Federal Systems Division. 

N \ The first study featured analysis of real satellite arrays, perform- 
ing spectrum analysis only, using the Fast Fourier Transform method. It 
was conducted by J. Leese, R. DePew, and B. Clark. The second study 
featured cross -correlation analysis of synthesized arrays using the lagged 
product method. This work was accomplished by B. Clark, J. Leese, and 
N. Woodrick. 

During the. remainder of 1968, work performed by B. Clark, C. Bevans, 
and J. Leese, resulted in the preparation of experimental computer programs 
to extract data in 64 x 64 arrays from ATS data tapes provided by the 
National Environmental Satellite Center and perform motion computations 
through cross-correlation analysis using the Fast Fourier Transform. Work 
performed by R. Stallard and P. Wrotenbery investigated the technique of 
applying Fourier Analysis and Covariance (or correlation) to cloud data. 

Since 1969, when he joined the organization, the application of the 
techniques to ATS data has been continued by Dr. J. A. Leese and his 
associates at the National Environmental Satellite Service. 
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The current study, which is the subject of this report, is a NASA- 
funded effort which investigates the feasibility of using these techniques 
to obtain cloud motion using data from a polar orbiting satellite. The 
Investigation has been performed by B. Clark of IBM with guidance and 
support from J. Greaves and W. Shenk of the NASA Goddard Space Flight Center. 

The major portion of the study was the development of computer programs 
to determine the amount of cloud motion during the interval between two 
overlapping satellite images. Debugging and initial testing of programs 
was done with array pairs previously generated from other source data. 

Early In the second quarter of 1971, upon instructions from the 
Technical Monitor, consideration was undertaken of the design of a compara- 
tive experiment in which the digital technique being Investigated under 
this contract would be compared to another method of estimating cloud 
motion by means of optical cross-correlation analysis, both using, ATS 
data. Further detail of this consideration appears in Appendix D. 

The investigation deferred agreement upon the type of data to be 
used (see Appendix D) for technique evaluation as well as the design and 
development of the array composition routines and acquisition of ground 
truth data. In the third quarter (September 1971) the Technical Monitor 
decided that the remainder of the study should be conducted as originally 
planned insofar as the remaining resources permitted. It was agreed that 
NASA would select data samples and provide data tapes to IBM within the 
succeeding two weeks. The number of tape pairs to be provided and the 
expected density of data points would have given from 15 to 50 cases with 
about 9 array pairs for each case, a total of 135 to 450 array pairs, avail- 
able by mid-October. 

Twenty cases were selected and requested by the Technical Monitor in 
the form of Grid Print maps, pseudo color photographs, black and white 
prints and digital tapes. Tape was the form needed for the testing of the 
programs in the study; the other standard output forms assisted in verifying 
and interpreting the tapes. 


1-2 



The first five sample cases proved unsatisfactory. From them it was 
apparent, however, that the source data did not include sufficient raw 
observations to yield useful data for a scale of 1 to 1 million. Hence 
further requests were for 1 to 2 million which would provide only about 
65 columns in the overlapped area. It developed that with this scale the 
maximum number of rows could be about 135 without excessive degradation 
except at two comers. 

When it became apparent that the remaining 15 cases could not be 
ready by the end of the eleventh month (January 1972) of the performance 
period a three-month contract extension was requested and granted. In the 
extended period four sample cases had been furnished by the start of the 
last month, two of which would form only one array pair each and two which 
would form five array pairs (not all independent but each pair displaced 
sixteen rows from its predecessor). These were used to test the array 
forming routines and the array processing programs; some program modifica- 
tions were made as a result of that testing. In the last month of this 
Extended period four more sample cases were furnished and the remainder 
of the twenty cases originally selected were found unsatisfactory. These 
last Vour cases were tried on several computer runs each but proved 
unreadable. A subsequent investigation revealed that the tapes were 
essentially blank due to the failure of the NASA programs used for the 
tape writing. 

A short extension of the performance period at slight added cost was 
then requested and granted. New versions of the last four sample cases, 
with extended latitudinal range, were received when work was resumed. The 
tapes were found readable and their data agreed with the data appearing 
on their corresponding grid print maps. 

In the limited time remaining available for computer processing some 
testing was conducted using the eight sample sets. Details describing these 
eight samples appear in Appendix E. A summary of tests performed appears 
in Appendix G while sample output results are included in Appendix H. 
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It became apparent as Is further described in Sections 5 and 6 that 
for THIR (and probably HRIR) data the applications of this technique 
to polar orbiting satellite data is impractical due to the difficulty of 
obtaining sample image pairs containing useful data. 
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Section 2 


SOME BASIC THEORY 

Two basic theoretical considerations in this study are the concept 
of cross correlation and the relation of cloud top temperature (or 
equivalent black body temperature) to cloud top height. 

2.1 The Cross Correlation Concept 


Two dimensional cross correlation provides the basis for the algorithm 
used to estimate motion which has occurred between two image arrays, 
observed at different times over the same geographic area. 

• 

Let x t denote the individual elements of a set X of measurements of 
some physical parameter and let the symbol E(T] mean the "expected value 
of Y". If the probability distribution f x (x) associated with the measure- 
ments is normal (Gaussian) it can be completely characterized by its mean 


V = E [X] 
and its variance 



2 2 /•” 2 
a = E [ ( X-/i) ] = / [x-fj) f y, (x) dx. 

— 00 

Then only for a purely random series will neighboring values be independent. 
In general, neighboring values of a time series will be correlated. Hence, 
in addition to specifying the mean and variance it is necessary in the case 
of a stationary normal series to specify its autocovariance function. For 
a specified lag A this may be expressed as: 


2-V 



c (X) = E [ (X(t) -U) • (X(t+A) -v) ] 

■ j (x t -/i) f x (x) dx 

- 00 

In practice C(X) can be estimated by 

N-X 

CU) = < x t+x '* ) 


where 


x = 


II 

1 £ X 

N x f 



When two processes, say W^ and Wg are being studied it is possible 
that they have different scales or different variances. In this c-ase 
we define the cross-correlation function for equal sized (say M-element) 
\jsets of elements from W-j and W ^ as 


R 12 (X) 


C (X) 

a l ° 2 


where C(X) 



and a-j and o 2 denote the standard deviations of the two processes. 


The essential features of two-dimensional cross-correlation are shown 
graphically in Figure 2-1. The input array consists of data values over 
an area for two different time periods a'2 depicted by X Q and X t in (a) of 
Figure 2.1. Cross-correlation coefficients are computed for different 
lag values of X^ upon X Q . The value at lag (p, q) is given by 


R (p, q) - 


Cov (p, q) 

O y CT y 

A o A t 


where R (p, q) Is the cross-correlation coefficient at lag values p and q 
in the P and Q directions, respectively. 
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c. Cross-Correlation Matrix 


Figure 2.1 Graphical Representation of Two-Dimensional 
Cross-Correlation 
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Cov (p, q) is the covariance at lags p and q, and 

°Xq and a X t are the standard deviations of the input arrays Xq 
and X t respectively. / 

The coefficients are determined over the limits 
-P < p< P 
-Q < q < Q 

as shown in (c) of Figure 2-1. The limits. of P and Q.are.a function of the 
particular problem and need not be symmetric in any direction. 


The locations of the maximum values in the cross-correlation matrix 
are related to the cloud motion in the following manner: 


where 


\ 



. [(p'Ax) 2 . ♦ (q-fly) 2 ] 

At 



p 1 AX 
q* Ay 


| C | is the speed of motion, 

6 is the direction of motion, 

p' and q' are the locations of the maximum in the correlation 


matrix 

Ax and Ay are the finite sampling Intervals for the input 
arrays in the x and y directions respectively. 

At is the time interval between the two pictures. 


Some supporting derivations may be found in Appendix I. 

■o . ’ 


o 
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2.2 The Relation of Cloud Top Temperature to Height 

The image data used in this study consists of radiometric observations 
derived from the 10.5 - 12.5 micron window channel of the two channel 
high resolution scanning radiometer used in the Nimbus IV Temperature- 
Humidity Infrared Radiometer Experiment. Based upon the sensor calibra- 
tion the radiance measurements have been converted to measures of equiva- 
lent black body temperature expressed in degrees K. The digital information 
is provided in this format, representing day or night cloud top or surface 
temperatures. These values of equivalent black body temperature have 
also been converted separately to units of grey scale which in turn have 
'been displayed as black and white photographic images. In some cases the 
temperatures have been expressed in a more detailed color coded form as 
colored photographic images. 

In using equivalent black body temperature data to represent -images 
in which clouds may be embedded, it must be recognized that: 

1. the temperature-height relationship is not biunique, 

2. there can be considerable overlap in temperatures which might 
be observed at cloud tops, snow, ice, land or water surfaces. 
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DATA FLOW 

The capability used in deriving cloud motion estimates is made up of 
several elements which were not all active concurrently during the course 
of this study. The description of the capability is therefore presented 
by the identification of the various elements followed by an indication, 
through the flow of data and of manual and machine aided functions, of 
how these elements contribute to the study procedure. 

3.1 Elements of the Cloud Motion Study Capability 

1. The Nimbus IV Satellite 

The Satellite was designed to be placed in an orbit which would 
be circular at 600 nautical miles, sun-synchronous having a local high 
noon equator crossing and an 81 degree retrograde inclination. Successive 
orbits would cross the equator at 26 degrees of longitude separation and 
the period would be about 107 minutes. It was launched successfully into 
a near circular orbit (at 587 x 593 nautical miles) on April 8, 1970. 

There were 10 meteorological experiments included in the spacecraft 
configuration. Details of the objectives and design characteristics of 
the Nimbus IV Spacecraft System appear in the Nimbus IV User's Guide. 

2. Temperature Humidity Infrared Radiometer (THIR) Subsystem 

The THIR Subsystem may be subdivided into three parts: 

(1) The on-board portion which in turn is made up of the two 
channel high resolution scanning radiometer, the High 
Data Rate Storage Subsystem (HDRSS), and the Real Time 
Transmission Subsystem (RTTS). The HDRSS and the RTTS 
are shared with other experiments. 
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(2) The ground based data collection portion where the THIR 
information is demultiplexed and recorded on magnetic 
tape. This is part of the Data Acquisition Facility (DAF). 

(3) The ground based data processing portion where the signal 
is demodulated, synchronized and recorded both as film 
strip and in digitized form on magnetic tape. This is 
part of the Nimbus Data Handling Facility (NDHF) located 
at Goddard Space Flight Center. 

Details concerning the THIR Subsystem are included in the Nimbus IV User's 
Guide and the volumes of the Nimbus IV Data Catalog. Some extracted 
information is also included in Appendix A. 

3. The NASA-Goddard Space Flight Center capability to prepare 
selected Nimbus IV sample data in the form of standard output products 
suitable for subsequent processing in the study. This involved both 
manual interpretation and selection procedures, the GSFC computing facility 
^nd certain computer programs. The sample selection was performed manually 
by\the contract Technical Monitor based upon film strips and related data. 
The Standard form reduced radiation data tape called the Nimbus Meteoro- 
logical Radiation Tape - THIR, abbreviated NMRT-THIR, was or had been 
generated on the IBM System/360 computer using routine procedures. The 
NMRT-THIR for each data sample was processed by NASA and contractor 
personnel using three programs. 

(1) The Nimbus HRIR Mapping program (NHM) which converts the 
digital data for a Mercator grid map. The Nimbus HRIR 
and THIR formats for the grid map are identical. The grid 
map is produced on a printer and the line images and related 
information are also output as unformatted FORTRAN records 
on a magnetic tape. 

(2) The Pseudo Color Input Tape Generation Program (PCITG) 
which generates a magnetic tape copy of the mercator data 
set placed on disk by the grid print mapping program NHM. 
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(3) The Pseudo Color Mercator Mapping Program (PCMM) which 
produces a magnetic tape in a format compatible with 
color photograph facsimile equipment to generate a pseudo- 
color mercator projection map of the THIR data. 

4. The cloud motion programs are also operated at the Goddard Space 
Flight Center computing facility under the sponsorship of the Technical 
Monitor. The first of these is an array formation program which reads 

data from two grid print output tapes, selects the desired records, organizes 
their data into subarrays and collates these for the two data times into 
array pair records on a single tape. The second is an array processing 
• program which performs cross correlation analysis and determines estimated 
vectors represented pattern motion. 

5. The Investigation function is performed by the Principal Investi- 
gator with guidance from the Technical Monitor. Results are compared to 
weather data obtained through the National Oceanographic and Atmospheric 

"x Administration (NOAA), estimated vectors are interpreted, and overall study 
results are evaluated by the investigator. 
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3.2 Flow Through the Elements of the Study 

The elements of the Cloud Motion Study procedure have been combined 
into a single diagram in Figure 3-1. The organization of the diagram 
presents the distribution of functions between these elements and suggests 
the routing of data through the major processing steps. Two diagrams from 
the Nimbus IV User's Guide provide more detail of the THIR Subsystem. 
Figure 3-2 shows the interrelation between the three portions of the sub- 
system in a simplified block diagram. Figure 3-3 further delineates the 
conversion from analog to digital data which is presented as a function 
of the Nimbus Data Handling Facility as well as the preparation of the 
NMRT-THIR tape which is the source of data for the Nintous HRIR Mapping 
program. 

In Figures 3-4 » 3-5 and 3_6 . are shown the broad functional flow 
through the Array Formation and Array Processing programs developed during 
the study. These two programs are combined with a control program in the 
.study to fit the polar orbiting type THIR data which was selected for the 
study. While the Array Processing program is generally adaptable to 
processing of pairs of rectangular arrays, the Array Formation Program and 
the Control Program which combine the two are dimensioned to fit the 
particular data used. 



3-4 





Figure 3-1. BROAD DATA AND FUNCTIONAL FLOW TOR STUDY 


3-5 



































ANALOG 



Figure 3- 3. Simplified Block Diagram of the A/D Processing System 

(From Nimbus IV User's Guide) 
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Section 4 


IMPLEMENTATION 

The overall approach to the study Included the selection of several 
previously developed computer programs, some incomplete, testing them 
with artificial image pairs, evaluating their various approaches and 
evolving therefrom a single procedure for deducing pattern motion. A 
type of experimental data was to be selected from polar orbiting satellite 
experiments and the necessary routines developed to prepare that data 
in format suitable for input to the motion computation procedure. .Actual 
data samples would then be used to evaluate the procedure and it would 
be modified accordingly. If resources permitted, statistical experimen- 
ts., tation would then be undertaken using a larger sample of data. 

4.1 Initial Technique Development 

Three experimental programs developed within IBM to compute pattern 
motion for meteorological satellite data had been used at three IBM 
System/360 Model 50 computer installations which used different compilers 
and operating systems. Other programs provided by Dr. J. A. Leese and 
Mr. C. S. Novak of the National Environmental Satellite Service (NESS) 
had been used at their Control Data Corporation system 6600 installation 
which, due to the different word length of the machine, had a slightly 
different language. All programs had been written in FORTRAN IV. The three 
IBM programs were adjusted to the Goddard Space Flight Center computer 
system and debugged. One of the NESS programs was converted to the IBM 
System/360 FORTRAN IV language and executed with data samples which were 
provided with them. 

The three IBM programs had been designed to serve distinct purposes: 
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1. The first was to test the Fast Fourier Transform method of cross 
correlation computation using small synthesized arrays. 

2. The second was a streamlined program which would perform 
only those functions necessary to derive a pattern motion 
vector and pHnt a minimum amount of information. 

3. The third was a more general program which would perform both 
cross-correlation and spectrum analysis and offer a more 
comprehensive printing capability. 

The first program was intended to test only the essential features 
which would be needed in the array processing program. Prewhitening of 
arrays before power spectrum computation, smoothing of the cross trans- 
form array before computing the cross -covariance, and filtering for 
selected wave numbers were not included as candidate functions. Inputs 
were constructed as 16 x 16 complex word arrays and a simple standard 
output format was used. Arrays were bounded with uniform strips equal to 
the array mean, wide enough to permit the assigned motion without loss 
of points whose value differed from the mean, and elements vacated during 
motion were reset equal to the mean. 

Features of the program were: 

a. Constructed ’one of five types of initial array and moved this 
with the specified vector to produce a second array. 

b. Transformed each array and computed a cross -transform array 
by multiplying each transform element of the second array 

by the complex conjugate of the transform element of the initial 
array. -o 

c. Formed the cross-covariance array by taking the Inverse trans- 
form of the cross -transform array. 

d. Computed the power spectra of both the input arrays. 

e. Located the highest cross-covariance value, found its lag 
coordinates and formed a relative cross -covariance array by 



dividing each cross -covariance by this maximum. 

f. Printed the array resulting from each of the nine array 
computations as well as the input type and motion vector. 

The second program was designed to perform the basic functions 
required to arrive at a motion vector from two input data arrays 
associated with the same geographic location’ but separated by a known 
interval of time; Features of the program were: 

a. Read in the data for each array and computed the mean and 
standard deviation of each then subracted the mean from 
each element of the respective array. 

b. Transformed these arrays and computed a cross -transform 
array by multiplying each transform element of the second 
array by the complex conjugate of the corresponding element 
In the first array. 

c. Smoothed the cross transform array with Hanning coefficients 
and formed the cross-covariance array by taking the inverse 
transform of the smoothed cross -transform array. 

d. Formed the cross -correlation array by dividing the cross- 
covariance array by the product of the two single-array 
standard deviations. 

e. Located the highest cross-correlation value and computed the 
motion vector corresponding to its location. 

f. Printed the time and identification for each input array, thei 
arithmetic means, the lag coordinates and value of the peak 
cross-covariance and the motion vector expressed in degrees 
and knots . 

The third program was designed to analyze each input array, compute 
a motion vector and print those results chosen from among the output 
options available, from a pair of arrays such as is processed by the 



second program. This program included all of the features of the 
second program and in addition, on option: 

a. Prewhitened each input array by making a least squares fit. 

b. Computed the power spectrum for each prewhitened array. 

c. Computed the cross-power spectrum of the pair of arrays. 

d. Printed those detailed results requested from these options. 

1) Print input arrays for both observation times in integer 
format. 

2) Print Fourier transform array for both observation 
times in complex format. 

3) Print power spectrum arrays for both observation times 
in integer format. 

4) Print the cross power spectrum relating both observation 
times in integer format. 

5) Print the smoothed cross covariance array In integer 
format. 

6) Print the cross\:orrelation coefficient as a percentage 
in integer format. 

The program developed at NESS processed 32 x 32 element arrays 
altering the array for the second time by keeping the central square, 
which contained one fourth of the input data for that time, and replacing 
the remainder of that array by the mean value of the array. It then 
performed the correlation analysis for the two resulting arrays. After 
this step it zeroed the portion of the transformed arrays corresponding 
to specified wave numbers, found the resulting filtered array and completed 
the analysis cycle for this new pair. The cycle could be repeated for 
up to five sets of wave numbers specified in input parameter cards. 



4-4 



4.2 Preliminary Timing Evaluation 

A preliminary timing evaluation was made using the second IBM 
program (referred to as the short program) and the third IBM program 
(long). The two programs were modified to allow the cross-correlation 
computation itself to be bypassed but to permit the remainder of the 
processing to be completed in the customary manner. Those functions 
which were bypassed were the performance of the Fourier transform of 
the input arrays, the cross product of one transform with the complex 
conjugate of the other, smoothing of the result using Hanning 
coefficients and taking the inverse Fourier transform of the smoothed 
array. The same set of input array pairs was then processed first by 
the original programs then by the modified by-pass programs to measure 
the performance in three ways: • 

1. Using the short program . 

2. Using the long program with -output limited 

3. Using the long program with full output 

The average times required for the computations, not including data 

s . 

selection and array formation ('since the input was preselected arrays 
in card format) were per pair of 64 x 64 element arrays: 


Program 

Option 

Cross - 
Correlation 

Input 

Output 

Other 
CPU Tima 

Total 

Time 

Short 

1.2 sec 

0.7 sec 

0.6 sec 

2.5 sec 

Long (limited) 1.2 sec 

1.2 sec 

1 .2 sec 

3.6 sec 

Long (full) 

1.2 sec 

2.1 sec 

3.9 sec 

7.2 sec 


This indicated that for the best case the cross correlation portion used 
about 48% of the total computation time for an array pair. 
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4.3 Technique Modifi cations 

1. Suppression of Spurious Products 

In computing a cross-covariance for two arrays of equal size 
for non-zero lags the sum of cross-products may be in error, if using 
the lagged-product method because there' are non-overlaoping portions, 
or if using the Fourier transform because of its implicit interpre- 
tation of both data arrays- as_repeating sets In each direction with 
period corresponding to their respective dimensions. Since this study 
) involves cross -correlation analysis using the Fast Fourier Transform 

this error is eliminated by using an array for the second time (search 
area) which is twice the desired size and modifying the array for the 
initial time (window area) by embedding it in a frame of zeros, after 
the mean has been subtracted from each element. The resulting cross- 
covariance array for lags corresponding to the interior of the zero- 
frame will have that error suppressed. 

'n. 

2. Restricting the Range of Data Values 

\ Experimental two-dimensional image data may be missing over 

certain data points or it may fiave values outside a reasonable range 
which would tend to distort cross -correlation computations. In the 
case of data from the Nimbus IV THIR Experiment, there is further risk 
possible because values which appear to be reasonable may have been 
measured over areas free of clouds. Missing data will be represented as 
zeroes and the lowest valid values are large; however, it is not always 

obvious that a value is too high. It was therefore necessary to provide 

means of screening out the effects of missing data and unwanted high 
data separately. This was done by accepting as program input parameters 
an upper and a lower threshold. 

In each array of a pair to be considered the individual data values 
which will enter into the cross-correlation computation are tested 
against the upper and lower bounds or thresholds. (In the case of the 

array for the initial picture only the central quarter or window area of 
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the array is involved whereas for the second picture the entire search area 
is considered.) Those which lie within bounds are used to compute a mean 
and standard deviation for the array. On a second search through each 
array, 

a. Data values which are too low are replaced by uniformly 
distributed random numbers whose mean is equal to and whose 
standard-deviation-isa' multiple of those of the corres- 
ponding array. This tends to eliminate the adverse effect 
of missing data. 

b. Data values which exceed the upper threshold are replaced 

by the mean value of the corresponding array. They thus 
do not contribute to the cross-covariance. . 

While the lower threshold value Is likely to be a function.of the 
sensor, the upper threshold may vary with season, time of day, or 
geographic location. In fact, the threshold approach offers the oppor- 
tunity to consider slicing even an image of valid information into distinct 
sublets for special analysis',. 

3. Normalization for Subsamples' 

Investigation of the subsample means and variances for each 
window-sized subarray of a search area array revealed significant varia- 
tion in values, sometimes by as much as a factor of 3 or 4 between the 
highest and lowest values of standard deviation. Since these standard 
deviations enter as divisors in computing the cross -correlation coefficient 
it is necessary to provide for effective normalization for each subsample 
and actual computation of the correlation coefficient prior to seeking 
the maximum value corresponding to the best pattern fit. A very 
efficient subroutine was programmed to accomplish this computation for 
all subsamples of the array for the second observation time. 

When the use of thresholds was coupled with the improved normalization 
for subsamples very dramatic effects occurred which complicated the 
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process of completing the cross -correlation computations and seeking 
their maximum value. In cases where there are very few clouds in the 
window area or in any search area subarray, the standard deviation 
computation may result in floating point underflow or division by zero. 
To remedy this, it was necessary to reject computations for such 
occurrences. The first fallback position in such a case is to use a 
constant subarray standard deviation equal to that of the search area; 
if that still fails computation is abandoned. 

4. Deletion of Non-essential Computations 


Several computer routines which had been considered in the 
program development were found non-essential to the cross -correlation 
analysis and were deleted from the program package in the interest of 
reducing storage requirements and computation time. Among these were 
the computation of power spectra for individual arrays or the cross 


spectrum for an array-pair, smoothing, which was useful in the spectrum 
computation but tends to dampen the extrema if applied to cross- 
correlation analysis, and the printing of Information not needed 
operationally. \ 


5. Simplification of Printing 


A new print routine was prepared to reduce printer output. 

The original output routine used four pages of printer output to display 
the array of cross-covariance or cross correlation values for the 
range of E-W lags from -32 to +31 increments and N-S lags from +32 
to -32 increments. The revised routine produces a 32 x 32 element page 
which is the window array extracted from the larger 64 x 64 array. It 
also prints the 32 x 32 subarray of the search area centered about the 
peak value of cross correlation. 
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4.4 Conversion of Input Sample Tape Pairs into Sets of Array Pairs" 

for Analysis 

Routines were programmed to read in each sample tape in a pair and 
subdivide its data into as many 64 x 64 element arrays as can be formed 
starting from the Northernmost row and Westernmost column and incrementing 
each starting row or column by 16 (a value which could be changed) before 
proceeding to the next pair. The .largest array. which may be provided 
on the sample tapes would be 200 x 200 data points; such a sample would 
yield at most 9x9 data arrays with three-fourths of the data elements 
being common in adjacent arrays. 


\ 

\ 
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4.5 Organization into Program Modules 

The main control program, CMXC, and the tape reading and array 
formation subprogram, RDTAPE, are necessarily written with dimensions 
which are specific to this study, computer system and type of data 
used as input. All of the remaining 20 subprograms have been prepared 
in modular form with variable dimensions to permit their adaptation 
to similar use with different data. 



Clock reading is introduced at the beginning and ending of each 
major operation and times used are printed out at each intermediate 
point. Summaries by each operator are also printed at the completion 
of the job. 



4.7 Characteristics of the Computer System Used in the Study 

The IBM System/360 Model 91 is an information-processing system 
designed for ultra high-speed, large-scale scientific and business 
applications. It provides a major-machine-cycle time of 60 nanoseconds. 
Data flow Is eight bytes (one double word) in parallel. The storage 
cycle time is 780 nanoseconds. (The cycle time of storage itself is 
750 nanoseconds.) Minimum total storage access time is 600 or 900 
nanoseconds, as determined by the way in which the processor storage 
is attached. 

Because floating-point overflow and underflow cause imprecise 
interruptions on the Model 91, It is possible that subsequent instructions 
will be executed using the overflow or underflow results. For this 
reason, the results are made to differ from the standard System/360 
results, which produce the correct fraction and a wraparound exponent. 

On the Model 91, overflow produces the correct sign and the maximum 
fraction and exponent; underflow produces a true zero result. 

The GSFC computer uses operating System/360 with MVT (Multi- 
programming with a variable number of tasks). Programs were compiled 
using FORTRAN IV, H Compiler, Option 2. 



4.8 Program Options 

To indicate the necessary parameters to the main program, there 
are nine input control card types. The options permitted are: 

1. Nunber of Input Tapes 

2 Indicates run is to process two NMR print tapes and generate 
an array tape before completing array processing. 

1 indicates an array tape Is used as the input tape. 

2. Nunber of Threshold Values 

The number of upper thresholds to be processed is entered. 

All computations wi 1 1 be repeated for each new one. 

3. Number of Replications 

For timing runs the number of repetitions of computations 
for each array pair may be specified. Normally it. is 1. 

4. Noise Scaling 

The number of sixths of array standard deviation to be used 
In computing random noise to replace missing data may be 
specified. Normally it is set at 6. 

5. Print Option 

1 indicates array printouts are requested 
0 Indicates they are bypassed. 

6. Initial Values 

The first upper and lower thresholds and an odd value for use 
in starting the random numben generator are specified. 

7. The run number and orbital description are entered to permit 
subsequent identification of results. 

8. Time between orbits (in minutes) if the run Includes array tape 
generation. 

9. Additional upper threshold cards if the number of threshold 
values to process exceeds 1. 



Section 5 


RESULTS 


5.1 Vector Computations 

The data sets used as input to the vector confutation programs are 
summarized in Figure 5-1. As the figure shows, eight sample pairs were 
provided as test cases. From these eight cases, a total of 32 array pairs 
were formed, as described in Appendix E. 

The 32 array pairs were processed using a value of 190 as the threshold 
for rejecting missing data. Three values (263, 268, and 273) of the 
threshold used to discriminate between cloud and ground target returns 
were used for each array pair. In addition, the five array pairs of 
Case IV-2 were processed against themselves as a control. 

Results of the vector computations may be classified according to 
whether a vector was apparently computed successfully; a warning appeared 
because the motion appeared to go to, or possibly beyond, the picture 
boundaries; or the computation was rejected because of the floating point 
underflow condition discussed on pages 4-7 and 4-8. Such a classification 
is shown in Figure 5-2. The results may be summarized as follows: 



Array Pairs 

Control Pairs 

Vector Found (V) 

53 

5 

Warning (W) 

7 

0 

Rejection (R) 

v36 

10 


As Figure 5-2 shows, there were a total of 21 array pairs (not 
Including control cases) for which at least one apparently valid vector 
was computed. For four of these array pairs the results were considered 
untrustworthy because less than 30% of the data values were considered 
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Figure 5-1. THIR DATA SAMPLES SELECTED FOR USE IN CLOUD MOTION STUDY 



















Case 

Array 
Pai r 

Upper Thresf 

old 

Comments 

263 

268 

273 

I 

1 

V 

V 

V 

. 

Data too sparse 

II 

1 

R 

R 

R 

i 


2 

R 

R 

R 



3 

R 

R 

R 



4 

R 

V 

V 

Data too sparse 


5 

R 

V 

V 

Data too sparse 

III 

. 1 

V • 

V 

V 



2 

V 

V 

V 



3 

W 

W 

V 



4 


V 

.V 



5 

R 

R 

R 

1 

IV 

1 

V 

V 

V 



2 

R 

R 

R 



3 

R 

R 

_ R 



4 

R 

R 

R 



5 

R 

R 

R 


V 

1 

V 

V 

V 

Insufficient ground truth 


2 

V 

V 

V 

Insufficient ground truth 


3 

R 

R 

R 



4 

R 

R 

R 

* 


5 

R 

R 

R 

* 

VI 

1 

V 

V 

V 

Data too sparse 


2 

V 

V 

W 



3 

V 

V 

V 



4 

V 

V 

V 



5 

V 

V 

V 


VII 

1 

V 

W 

W 


VIII 

1 

R 

w 

V 



2 

V ' 

V 

V 



3 

V 

V 

V 



4 

V 

V 

V 



5 

V 

V 

W 


Control 

1 

V 

V 

V 



2 

R 

V 

V 



3 

R 

R 

R 



4 

R 

R 

R 

... 


5 

R 

R 

Tt 



V = Vector Found W = Warning R = Rejection 


Figure 5-2. CLASSIFICATION OF VECTOR COMPUTATION RESULTS 


5-3 





























-£S63? f 

valid. For two additional array pairs available' ground truth was insuffi- 
cent to permit evaluation of the results. The |five vectors comouted in 
the control trials ^ showed no motion (as expected) ,< and all but one of the 
rejections occurred with less than 10 percentof^the data values acceptable. 

The computed results for the remaining fifteen array pairs were 
compared to rawinsonde data for the nearest observation time which could 
be found in the archives of the Atmospheric Sciences Library of the National 
Oceanographic and Atmospheric Administration. Observing stations within 
the image area for each case were identified and data obtained for each 
station. Those stations lying within the array boundaries were examined 
and based upon the approximate elevation in that area of the center of 
coldest equivalent black body temperature found in the image a representa- 
tive wind observation was selected. These are presented in Figure 5-3. 
Because of the breadth of the image area, and the large difference, between 
the times of the image sensing and the rawinsonde observation. It is 
"^difficult to evaluate this comparison; it is significant that in several 
examples there appeared fairly close agreement. In most examples visual 
inspection of the image photographs substantiated the computed vectors; 
in several cases the photographs revealed strongly circulating patterns 
which would have both rotational and translational components . 

In the sample display of test results shown in Appendix G are included 
four samples of results obtained from Case VIII. Figure G-2 shows that for 
the subarray centered at 45.0 N latitude and 65.0 W longitude, over 
Newfoundland, the window area had a standard deviation of 0.577. For one 
or more subarrays of the search area the subarray standard deviation was 
either zero or so small that its product with Jthat of the window area 
resulted in a value less than 0.01. The sub array standard deviation was 
then reset equal to the search area standard deviation and the problem 
still existed. Apparently for the upper threshold of 263 there are so few 
data that depart significantly from their me ariX value that no valid cross 
correlation array can be determined. On the other hand, for the same 
sample pair the array pair centered at 50.0 N latitude and 65.0 W longitude 




Figure 5-3. RESULTS OBTAINED FOR THE FIFTEEN SELECTED ARRAY PAIRS 
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determines the same vector for the three thresholds tested. In those cases 
the window subarray contains over 97 % accepted data points for the three 
tests. The window area appears to cover the southern and colder half of 
a cyclonic circulation which sweeps across the search area from West to 
East. These examples show the great differences in computational results 
that can be obtained from different array pairs drawn from the same sample 
pair. 



5.2 Timing Results 

Timing was estimated by. making a run using a five array-pair 
sample with the same four threshold values and replicating the compu- 
tations 20 times each. This provided 400 vector computations, many of 
which would find subarray standard deviations too small and result in 
recycling to attempt to find a vector assuming a constant scaling in 
place of the subarray standard deviations. The largest variation In - 
the operations was found In operation 1. Tape reading and array ... 
formation need be done only once for a picture pair so Its timing was 
measured independently. 


The typical times required 
Tape reading and array formation 


for a 1 -array pair 

0.3 

to 0.5 sec 

for a 5-array pair 

0.4 

to 0.7 sec 

Operations per array pair 



1. Initialize to process one pair 


0.58 sec 

2. Fourier transform array pair and form 
cross-product 


0.27 sec 

3. Inverse Fourier transform 


0.21 sec 

4. Rearrange cross covariance and find 
motion vector 


0.02 sec 

5. Printing (only If requested) 


0.40 sec 

Pair total with printout 


1.48 sec 

Pair total without printout 


1.08 sec 


The Initializing to process one pair which includes the threshold 
testing and subarray standard deviation computation uses about 54% 
of the time required to process a single array pair. 
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5.3 System Versus Cloud Motion 

The geographic scale of data points used In the study leads to 
a situation which had not been anticipated originally. At leas three 
of the picture pairs, making up 11 of the 16 useful arrays, though 
chosen because they were .in overlapped portions of successive orbits 
and between selected latitudes because they would otherwise have been 
severely truncated on the corners, were found to have apparent closed 
circulation patterns lying within their central portions. One instance, 
Case VII, proved to be a remarkable case which occasioned the issuance 
of two of the warnings because the cross correlation peak appeared at 
a boundary. 

Case VII is a single array pair case located over Southeast Siberia. 
The motion computations, one of which is presented as Figure 5-4, showed 
vectors from 7 degress with speeds of 41, 43 and 43 mps . The motion 
was measured between pictures completed at about 01:14 and 03:02 
^Greenwich Mean Time. 

Ground truth 400 mb level data taken at 000:00 Greenwich Mean Time 
showed a wind of 260 degrees at 55 mps at the center of the picture 


However, other winds in the picture area were: 

NW Comer 

360/65 mps 

N Center 

045/22 mps 

NE Comer 

195/24 mps 

SE Comer 

255/15 mps 

S Center 

285/23 mps 

SW Center 

295/56 mps 

W Center 

325/47 mps ^ 


The ground truth temperature, humidity, height and wind data for the 
400 mb constant pressure surface was plotted and contours of constant 
height drown. This is presented as Figure 5-5 where the computed motion 
vector appears as a dashed vector of 42 mps. Apparently intense northerly 
flow swept into the picture by the time the satellite passed overhead in 
its second orbit. 
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This case is a good example of the difficulty encountered In comparing 
computed vectors with ground truth data. Significant changes in the 
direction and intensity of motion over an image area can occur without being 
detected in the time between successive satellite orbits; to be aware of 
these it is necessary to study a synoptic analysis of the area and its 
environs. Motion such as was computed for this case could be valuable in 
making upper air analysis but there is a risk in an automated procedure 
that such results will be rejected as inconsistent and the information thus 
lost. The motion detected in this instance revealed short term changes 
which could be missed in the conventional observation system. 




5-10 



Figure 5-5. HEIGHT AND WINDS (IN WS) AT 400 rt SURFACE AT 0000 GMT, 
17 MAY 1970 - CASE VII. 
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Section 6 


CONCLUSIONS AND RECOMMENDATIONS 

6.1 General 

The primary conclusion to be drawn from this study is that the basic 
technique of cross-correlation using the Fast Fourier Transform can be 
applied to data from a sensor such as the THIR to determine cloud motions. 
Unfortunately, the small number of cases for which both computed wind 
vectors and ground truth measurements were obtained makes this conclusion 
somewhat tentative. It is felt, however, that the results presented in 
Section 5 do support the validity of the. process. 

\ The great difficulty in obtaining useful data leads to another, 
perhaps more important, conclusion: the combination of a polar orbiting 

satellite and a sensor such as the THIR is not practical for wind velocity 
determination. The sensor characteristics, limited sidelap of successive 
orbits, and present NASA procedures for generating geographically registered 
data sets make it quite difficult to obtain a sufficient sample population 
to support the required computations. Excessive manual involvement in 
sample selection and the number of steps needed to obtain digital output 
clearly show that it would not be economical to make cloud motion computa- 
tions with this type of data on an operational basis. 

A problem, common to all cloud motion determination techniques, which 
was observed during this study is that of discriminating between data 
samples that represent clouds and those that arise from ground targets 
such as ice or water. Simple thresholding was used in this study to 
accomplish the required discrimination. Multispectral signature analysis 
and change detection are examples of alternate techniques which might be 
used for this function. 
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6.2 Program Changes 

During the course of the study, the following changes were found to 
greatly improve the usefulness of the program in processing THIR data: 

a. Suppression of spurious cross covariance products by bounding 
the window area with neutral values (the window mean). In 

the experimental programs window and search area sizes were equal, 
whidrhad the effect of moving false clouds over the search area. 
The revised version measures the motion of only those clouds 
which lie within the window area initially. 

b. Replacement of missing data with noise on the basis of a lower 
threshold test. Missing data appears on the sample tapes as 
zero values. Since the lowest value expected in the samples 
would be about 190, missing data would have biased results 
toward the low side. By replacing missing data with random 
noise within the expected range this effect is neutralized. 

s..., 

c. Normalization of cross correlations according to the subsample 
scale (standard deviation). The variation in subsample standard 
deviation is very large. The original program found the cross 
covariance which is the cross correlation before it has been 
normalized to fit the scale of the two subsamples (the window 
and the search area subarray for that lag pair). Searching for 
the maximum of cross covariance would yield a good fit only if 
the search area subarrays had about the same scale. 

d. Exclusion, using an upper threshold test, of data which while 
valid to the sensor constitutes^noise to the cloud motion 
computation. The most difficult aspect of objective interpre- 
tation of the THIR data is in determining what is valid data 
for the interpretation. It is possible to receive measurements 
throughout the range of the sensor, but for motion study there 
must be some way of determining what is moving and what is fixed. 
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Temperature may be a good Indicator, though its effectiveness 
does depend upon sensor and geographic location. The uooer 
threshold provides a means of screening out data which is 
obviously too high. A summary of the effects of using three 
different threshold values appears in Figure 6-1. 



Upper Threshold U^per Threshold Upper Threshold 

= 273 / =268 = 263 


1 


2 


3 


4 


5 


Array Pair 


Case I 

2 

Single Array Pair 



II 

0 

0 

* 

3 

13 

III 

27 

31 

17 

16 

5 

. IV 

57 

28 

4 

0 

0 

V 

32 

42 

27 

* 

0 

VI 

5 

20 

64 

94 

82 

VII 

43 

Single Array Pair 



VIII 

16 

64 

98 

99 

96 

Control 

68 

38 

5 

0 

0 

Case I 

2 

Single Array Pair 



II 

0 

1 

4 

8 

21 . 

III 

33 

40 

31 

30 

11 

IV 

62 

35 

8 

0 

. 'o 

V 

38 

49 

30 

* 

0 

VI 

7 

25 

71 

99 

88 

VII 

51 

Single Array Pair 



VIII 

20 

v 70 

t 

v 99 

99 

99 

Control 

76 

43 

8 

0 

0 

Case I 

3 

Single Array Pair 



11 

* 

3 

8 

14 

28 

III 

43 

53 

45 

19 

25 

IV 

66 

40 

11 

0 

* 

V 

46 

56 

32 

1 

0 

VI 

16 

36 

78 

100 

94 

VII 

59 

Single Array Pair 



VIII 

30 

77 

99 

100 

100 

Control 

81 

48 

10 

0 

* 


Figure 6-1. DENSITY IN PERCENT OF WINDOW ARRAY POINTS QUALIFYING AS 
CLOUDS FOR THRESHOLDS SHOWN (* INDICATES LESS THAN U 
BUT NOT 0) 
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It appears reasonable to compute motion vectors using a 32 x 32 
window over a 64 x 64 search area. Time required for one such comouta- 
tion is about 1.1 seconds. 


6.4 Recommendations 

The obvious recommendation arising from this study Is to apply cloud 
motion determination techniques to data gathered by geosynchronous 
satellites. Use of such data should greatly reduce, if not eliminate, the 
problems of obtaining timely ground truth measurements and sufficient 
data sample population in overlapped areas. Data from both the visible and 
infra-red ranges of the spectrum should be used in such an investigation. 

Another area recommended for investigation is that of separating 
returns from clouds and ground targets. A simple approach that might be 
tried is that of differencing successive data sets and filtering out all 
'areas for which no appreciable change has occurred. An alternate approach 
is that of examining spectral and thermal signatures in order to arrive 
at an algorithm which will accurately detect clouds and reject ground 
targets in a single data set. 

A final recommendation is that a cloud motion determination scheme 
based on the Sequential Similarity Detection Algorithm (SSDA) 
be investigated. Recent experiments have shown that, especially for low- 
resolution data of the type used in cloud motion studies, use of the 
SSDA may produce results which are computationally much less expensive 
than those produced with cross correlation and the FFT with no loss In 
accuracy. 


1. Bamea, D.I., and Silverman, H.F. , "The Class of Sequential Similarity 
Detection Algorithms (SSDA's) for Fast Digital Image Registration," 

IBM Research Report RC-3356, May 10, 1971. 

2. Bernstein, R., and Silverman, H., "Digital Techniques for Earth Resource 
Image Data Processing," IBM Report FSC 71-6017, September 30, 1971. 

'O 

3. Barnea, D.I., and Silverman, H.F. , "A Class of Algorithms for Fast 
Digital Image Registration," IEEE Transactions on Computers, Feb. 1972. 
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Section 7 


NEW TECHNOLOGY 

No reportable items, as defined in Section 1(a) (i) of NASA Form 1162, 
NEW TECHNOLOGY (May 1966), have been identified during the performance 
under this contract. . 



Appendix A 
SENSOR DESCRIPTION 

(Condensed From Nimbus IV User's Guide) 

The THIR is a two channel high resolution scanning radiometer designed 
to perform two major functions. First, a 10. 5-12. 5/x window channel provides 
both day and night cloud top or surface temperatures. Second, a water 
vapor channel at 6 . 7/x gives information on the moisture content of the 
upper troposphere and stratosphere and the location of jet streams and 
frontal systems. The ground resolution at the subpoint is 8 Km for the 
window channel and 22 Km for the water vapor channel. The window channel 
"■'\will operate day and night while the water vapor channel will operate 
mostly at night. (Data from the window channel was used in this study.) 

The THIR radiometer consists of an optical scanner and an electronic 
module. A blackened collar near the scan mi r row serves as a sun shield 
to prevent sunlight contamination during spacecraft sunrise and sunset. 

The other side of the sun shield Is painted white. The end of the scanner 
opposite the sun shield contains the optical system and preamplifiers. 

The optical system consists of a scan mirrow, a telescope (comprised 
of primary and secondary mirrors) and a dichroic beamsplitter. The scan 
mirrow, inclined at 45° to the optical axis, rotates at 48 rpm and scans 
in a plane perpendicular to the direction^ of the motion of the satellite. 

The scan mirror rotation is such that, when combined with the velocity 
vector of the satellite, a right-hand spiral results. Therefore, the field 
of view scans across the earth from east to west In daytime and west to 
east at night, when the satellite is traveling northward and southward 
respectively. 
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The telescope focuses the energy at the dichroic beamsplitter which 
divides the energy spectrally and spatially into two (2) channels. A 
21 -mi Hi radian channel detects energy in the 6.7 micron band. A 7.0 
milliradian channel detects energy in the 10.5-12.5 micron band. It 
consists of a bandpass filter (transmission portion of the dichroic), an 
Itran-2 relay lens which also serves as a long wavelength blocking filter, 
a folding mirror, and a germanium immersed thermistor bolometer. 

The signals from the detectors are capacitively coupled to the pre- 
amplifiers, amplified and sent to the electronic module. In the electronic 
module, the signals are further amplified and corrected for detector time 
constant to provide the overall frequency response as required by the 
subsystem optical resolution. The signals are processed out of the 
electronic module through buffer amplifiers. The 6.7 micron channel -output 

is available on a full time basis as the shifted level channel. The offset 

* 

of the shifted level channel is provided in the buffer of that channel. A 
second video output selects either the 6.7 micron or the 11.5 micron 
channel by means of a command relay. In addition to the two video output 
signals, there are fourteen telemetry channels: ten analog and four digital. 

The instantaneous field of view (IFOV) of the window channel is approxi- 
mately 7 milliradians. At an altitude of 1112 kilometers (600 nautical miles) 
this results in a subsatellite ground resolution of 6.67 kilometers (4.1 
nautical miles). The scan rate of 48 rpm produces contiguous coverage along 
the subsatellite track. Due to the earth-scan geometry of the THIR, as nadir 
angle increases, overlapping occurs between consecutive scans, reaching 350 
percent overlap at the horizons, and resulting in a loss of ground resolu- 
tion in the direction of the satellite motion. Even greater loss of resolu- 
tion occurs along the scan line (perpendicular to the line of motion of the 
satellite) because of the expansion, with increasing nadir angle, of the 
target area viewed. 

Figure a_i shows, for the window channel, the relationship between nadir 
angle and ground resolution element size along the path of the satellite and 
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(b) 

Figure A~l. Relationship between Nadir Angle and Ground Resolution for the THIR 11.5/x. 
Chonnei at 600 N. Miles (o) Pictorial (b) Graphical 




perpendicular to it. The numbers under each resolution element are nadir 
angle (in degrees), resolution along the scan line (in kilometers), and 
resolution parallel to the satellite line of motion (in kilometers). 

No image is formed within the radiometer: the THIR sensor merely 

transforms the received radiation into an electrical (voltage) output with 
an information bandwidth of 0.5 to 360 Hz for the 10.5-12.5 micron channel 
and 0.5 to 120 Hz for the 6.7 micron channel. The radiometer scan mirror 
continuously rotates the field of view of the detector through 360 degrees 
in a plane normal to the spacecraft velocity vector. The detector views 
in sequence the in-flight black body calibration target (which is a part 
of the radiometer housing), outer space. Earth, outer space, and returns 
again to intercept the calibration target. The space and housing- viewed 
parts of the scan, which can be identified without difficulty, serve as 
part of the in-flight check of calibration. Information on housing tempera 
ture, which is monitored by thermistors, is telemetered to the ground 
stations and for calibration purposes is constantly compared with the 
temperature obtained from the radiometer housing scan. Even though the 
first stages of amplification are capacitor-coupled, the low frequency 
cutoff is so low that a dc restore circuit Is necessary to provide a zero 
signal reference. This occurs during that portion of the scan when the 
optics are receiving zero radiation (space). The dc restore circuitry also 
provides additional gain to raise the signal to the desired output level 
and filtering to establish proper frequency characteristics. 
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Appendix B 


DATA GENERATION PROCEDURE 


Following are summary descriptions of the three programs used at the 
Goddard Space Flight Center computer facility to generate the standard 
sample products provided for this study from the Nimbus Meteorological 
Radiation Tape. 

1. Program Number - L00240 (NHM) 

Title - Nimbus HRIR Mapping Program 

Abstract - , 

This program is used to generate maps of the Earth showing high 
resolution infrared radiation measurements taken by the NIMBUS NMRT-HRIR 
scanning equipment. Up to three maps can be made during a single pass 
of the NMRT-HRIR tape: - one mercator map and two polar maps. Thus, the 

entire Earth or any portion of It can be mapped at one time. 

Restrictions - 

Maps are limited to a width of 100 horizontal grids (25 grids 
per page). There is no limit to the number of vertical grids. 

The only limit to the overall map size is available memory space. 
Four bytes are required per map grid for each of the maps. 

2. Program Number - S00009 (PCITG) 

Title - Pseudo Color Input Tape Generation Program 

Abstract - 


This program is to be Included as an additional job step following 
a job step using program L00240. The HRIR (and THIR) grid print mapping 
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program (L00240) generates a two dimensional matrix containing tempera- 
ture data geographically located on a mercator map projection and places 
it (the matrix) on disc for passage to another step (L00244) which formats 
the matrix for printer (Class M or A) output. S00009 accesses the matrix 
on disc and copies it with necessary documentation, on magnetic tape 
(using an unformatted write). 

Restrictions - 

Limitations are set by L00240, and S00010. 

Current limit is 200 by 200 for Data and Population Matrices. 

3. Program Number - S00010 (PCMM) 

Title - Pseudo Color Mercator Mapping Program 

Abstract - * 

S00010 is a program for producing pseudo-color maps of data from 
the NIMBUS 3 HRIR and the NIMBUS 4 THIR experiments (or any compatible 
"'data source). The program accepts a data matrix from magnetic tape (see 
S00009) and reformats the output magnetic tape line by line to produce 
the pSeudo-color map containing a^ title, an annotated color scale, geo- 
graphical gridding, and the pseudo-color data map in a mercator projection. 
The gridding and associated annotations may be omitted if a non-mercator 
projection data matrix is used. A variable number of colors (up to 20) 
are available in the color scale. 

Restrictions - 

Data Array maximum size (MRCAT Tape) is 200 x 200 (1*4) Words. 

Data Values are limited to range of 1 to 350. Values ^ o are not contoured 
and represented as Black on Color Picture (output). Values > 350 are 
revalued to 0. Maximum of 20 colors can be selected. 
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Appendix C 

COMPUTER SYSTEM DESCRIPTION 

The portion of the Goddard Space Flight Center computing facility 
which was used for the cloud motion program development and testing during 
the study was the IBM System/360 Model 91K located in Building I. Support- 
ing services at the computing facility were also available when needed. 

The system hardware, on-line and off-line peripheral devices and supporting 
services are described in the SESD Computing Center User's Guide, includ- 
ing a diagram of the system configuration. 

The main system components include: • ’ 

a. Model 2091K Processing Unit, including one Model 2250-1 Display 
Unit serving as the operator's console. 

b. One Model 2395-1 Processor Storage, 2048K (2,097,152) bytes of 
high-speed storage. 

c. One 2150/1052 typewriter console. 

d. Three 2860 selector channels with: 

1. Two 2314 Direct Access Storage Facilities, containing 
233.408K bytes each. 

2. Two 2301 Drum Storage, containing 4000K bytes. 

3. One 2321 Data Cell, containing 400.000K bytes. 

e. One 2870-1 multiplexer channel, including three selector sub- 
channels with: 

1. One 2250-1 Display Unit 


2. Six 2401-3 7-track tape drives 

3. Eight 2401-6 9- track tape drives 

4. Two 2540 Card Read Punch 

5. Four 1403 N-l Printers 

6. One 2702-1 remote comnuni cations device for attaching 
Model r 1050 CRBE terminals. 
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Appendix D 


CONSIDERATION OF AN ALTERNATE EXPERIMENT 

Early in the second quarter of the study (the summer of 1971) upon 
instructions from the Technical Monitor, consideration was undertaken 
of the design of a comparative experiment in which the digital technique 
being investigated under this contract would be compared to another 
method of estimating cloud motion by means of cross -correlation analysis, 
both using ATS data. 

The instructions followed in this consideration initially were to 
consider the design of a comparative experiment, using ATS data, in which 
two methods of estimating cloud motion by means of cross -correlation 
analysis would be compared. One technique was to be optical/analog and 
the other the digital technique being investigated by IBM under Contract 
NAS5-11859. Subsequently these instructions were changed; instead of a 
comparative experiment, a comparative test would be considered and the work 
would include recommending a suitable set of test ATS data and reorienting 
the conduct of the contract performance to the use of ATS data instead of 
polar-orbiting data. 

Investigation suggested that the most serious problem in such a test 
might be that of registration relative to a frame of reference. It seemed 
unlikely that registration for ATS pictures will be automated in the near 
future. Manual performance of the task at the National Environmental 
Satellite Service had Indicated that the validity of the result diminishes 
with distance from the landmarks used for reference. Solution of that 
problem was considered well beyond the scope of contract NAS5-11859. 

A second problem would be that of data selection, availability and format. 
ATS meteorological data processing was being conducted by NESS so NASA did 
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not maintain complete data catalogs. 

To address these two problem areas it appeared possible that test 
data be selected from that previously . registered and processed by the 
National Environmental Satellite Service for the GARP month of June, 1970. 

In this way it would have been possible to eliminate the registration 
portion of the task which seemed likely to remain largely manual in the 
immediate .future, Independent of the method used to perform cross- 
correlation analysis. 

It was concluded that while consideration of the conduct of a compara- 
tive experiment or test using ATS data and the reorientation of the 
contract performance toward the use of ATS data as input would have been 
consistent with the broad aspects of the Statement of Work, it would have 
rendered some details of that statement inappropriate and would have called 
for a substantial departure from the proposed approach to this investigation. 

It was therefore recommended that the cloud motion study be completed 
as originally planned using HRIR data obtained from polar-orbiting NIMBUS 
satellites, that any comparative experiment or test conducted under the 
existing contract be based upon such HRIR data and that adaptation of 
programs and techniques to ATS data be considered separately. 
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Appendix E 

DATA SAMPLES SELECTED BY NASA 

Selection of data samples to be used in the study was accomplished 
by the Technical Monitor. The objectives of the selection process were 
that the sample pairs selected would: 

number between 15 and 50 

include about 100 by 100 data points each 

be selected from higher latitudes from overlapping frames 

of successive NIMBUS orbits of either HRIR or THIR data. 

The selection process included investigation of the NIMBUS IV Data 
^Cqtalog and examination of enlarged photographic prints of the Tempera- 
ture^H umi di ty Infrared Radiometer montages. Selections were made from 
NIMBUS"' IV 11.5/4 THIR daytime data. 

As pointed out in Volume 2 of the NIMBUS IV Data Catalog, the 
quality of THIR data from the window channel (11.5 /cm) was excellent. 
Data recorded after orbit 450 did contain periodic noise which could be 
filtered out. However, difficulties in processing to obtain the 
mercator tapes, sparsity of data which forced use of a coarse scaling 
and failure to find satisfactory data from both orbits chosen in the 
overlapped area which had been demarked resulted in degradation of the 
twenty sample sets (from forty orbits) into only eight useful sample 
pairs. These are identified as. Cases I through VIII which are described 
in Figure E-l. 

Three samples have been pictured to identify some of the limitations 
encountered in the available sample data. Figure E-2 shows the data 
selected for Case VI. As will be noted, a substantial field of cyclonic 
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circulation is nearly centered in the overlapped area; the limit of the 
East-West range to one array width restricts the analysis of variations 
in the field to those in the North-South direction. A central array pair 
is likely to provide information about motion of the field itself rather 
than within the field 

Figure E-3 portrays the data selected for Case I. Examination of 
the data values contained in this overlapped area reveal that only a very 
small portion of the central quarter of the area contains any cloud cover 
at all. This case was limited to a single array in both latitudinal and 
longitudinal directions so no significant conclusions could be made for 
this case. 

Figure E-4 portrays the data selected for Case VIII. In this case, 

the cloud features seem fairly well distributed with less prominent 

circulation features than were found in Case VI. Here the lower left 

♦ 

corner of the image from orbit 702 and portions of the right side of the 
image from orbit 703 are characterized by missing data values. 
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Figure E-l. THIR DATA SAMPLES SELECTED FOR USE IN CLOUD MOTION STUDY 
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!. Black and White Presentations of THIR Data Selected for Case VI 










NIMBUS 4 11.5 MICRON CHANNEL DATA 
D/Q 0403 05/08/70 08:02-08:07 GMT 1S2H FOR J. GREAVES 


337- W 305*8 303*W 30I*W 2WK 297*8 293*8 293*W 



.III ..I. .. .1.1 .1.1 .1 .1 1 


313 K 
300 K 
299 K 


29* K 
290 K 


284 K 
280 X 
279 K 

273 K 

274 K 
270 K 
269 K 

263 K 

264 K 
260 K 

239 K 

233 K 

234 K 
230 K 
249 X 

243 K 

244 K 

240 K 

239 *C 

233 K 

234 K 

° 230 K 

229 K 

223 K 

224 K 
220 K , 

; 219 K | 
• 213 K 

214 K 
' 210 K 
209 K 
203 K 1 



Figure E-3. Black and White Presentations of THIR Data Selected for Case I 





Figure E-4. Black and White Presentations of THIR Data Selected for Case VIII 



Appendix F 


COMPUTER PROGRAM 

The MAIM program which controls the flow of processing in the cloud 
motion cross correlation problem has been assigned the identifying label 
CMXC for use in referring to it as a subprogram. On subsequent pages 
of this appendix are the detailed source code (in FORTRAN IV Language) of 
the 22 subprograms developed in this study. Activity is actually divided 
into six basic operations each of whose central processor usage is timed 
by CMXC using the internal clock. Several subprograms are used in more 
than one operation. 

1. Tape Reading and Array Formation Operation 

The first time a pair of image tapes is processed an array pair 
tape is generated for use in subsequent operations using the subprograms 
RDTAPE, RDREC, RDPIK and RAYSET. For all comDuter runs the subprogram 
HARM is here initialized to process the appropriate size of arrays using 
subprograms SETHRM, HARM and IFEXIT. 

2. Operation 1 

Working storage is initialized to process one pair- of arrays and 
each of these arrays undergoes single arr^ analysis and, if necessary, 
modification to assure its data content meets the boundary conditions set 
for this pass through the loop. The subprograms used are NULOOP, SETRAY, 

Z FRAME, ZNORM, TIMNSD, THRESH and BKGRND. 

3. Operation 2 

The two modified input arrays undergo a Fourier transform and 
the resulting transformed arrays are combined into a product array W 
using subprograms HARM, IFEXIT and SPLITV. 
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4. Operation 3 

The inverse Fourier transform of array W is taken to determine the 
cross covariance array for the pair of arrays using subprograms HARM 
and IFEXIT. 

5. Operation 4 * 

The cross covariance array is rearranged in lag form and the 
corresponding array of cross correlation coefficients is computed. This 
array is searched for a peak positive value which, if found, is then 
converted into the motion vector which it signifies. Subprograms LAGWS, 
WINDOW, XCMAX, XIJMAX and VECTOR are used in this operation. 

6. Operation 5 

If requested in the run parameters, selected arrays of input and 
computed data are printed using subprogram OUTLAG. 



MAIN PROGRAM FOR COMPUTATION OF CLOUD MOTION 
COMMON INTEGRI 88192 ) 

DIMENSION U< 64,64 ) ,VC64, 64) ,W( 64,64 ) ,Z (64,64) ,WT I 65, 65 ) 

A ( 64,64 ) ,B( 64,64) ,C( 64,64 ) , TARA ( 64,64) , T ARBI 64,64) 

YWMNI 33,33), YWSD ( 33 , 33 ) , XCVP ( 33 , 33 ) , VARA ( 3 2, 32 ) 
XCV(33,33) ,S(64) ,RL( 5) 

MA(64,64),MB(64,64),KARA(t>4»64),KARB(64,64)»LARA(32t32) 
INV(64),LIST(32) »MH( 3) 

L AGRA Y I 33, 33 ) , AXCLDI 18 ) , TI M( 5, 500 ) 

REPL I 18 ) ,THR( 18 ) ,TAPEN( 18 ) ,PRNTN< 18) tSDEVNI 18 ) ,RN( 18) 

U, V,W,Z,WT,ZUVW 
ZI 

1), INTEGRI D) 

1 ) , INTEGRI 16385 ) ) 

1 ) , INTEGRI 32769) ) 

1 ) , INTEGRI 45315 ) ) 

1 ) , INTEGRI 53507 ) ) 

1 ) , INTEGRI 61699 ) ) 

I XCVP I 1, 1 ), INTEGRI 63877 ) ), I VARA I 
(XCV I 1,1),INTEGR(65990) ),(TIM I 
I SID , INTEGRI69579) ) 

IMA I 1, 1 ), INTEGRI 69648 ) ) 

I KARA! 1,1), INTEGRI 77840) ) 

I LARA I 1,1), INTEGRI 86032 ) ) 

I L I ST I 1 ) , INTEGRI 87120) ) 

(7X,I4,9X,I4,26X,I 10) 

( 10H0LATITUDE=»F9.5» 13 HN 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
COMPLEX 
COMPLEX 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
\ EQUIVALENCE 

"EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

1 FORMAT 

2 FORMAT 


I U ( 1 , 

I W( 1 , 

I WT ( 1 , 
( 8(1, 
( TAR A I 1, 
( YWMNI 1 , 


( VI 
( Z( 
( A ( 
( C( 
( TARBI 
( YWSDI 


I RLI 
( MB ( 
(KARBI 
( 1 NV ( 


1, 

1* 

It 

It 

It 

It 

1 

1 

1 

1 

1 

1 


1 ) 
1 ) 
1 ) 
1 ) 
1 ) 
1 ) 


INTEGRI 8193)) 
I NT EGR I 24 5 77 ) ) 
I NTEGR (41219) ) 
INTEGRI 4941 1 ) ) 
INTEGRI 57603) ) 
INTEGRI 62788) ) 
INTEGRI 64966) ) 
I NTEGR I 67079 ) ) 
I NTEGR I 69643 ) ) 
INTEGRI 73744) ) 
INTEGRI 81936) ) 
INTEGRI 87056) ) 


1ZE=,F8.5,19HDEG N-S 
2F6.0 1 8HMINUTES ) 

3 FORMAT 1 1 H 1 ) 


L0NGITUDE=,F9.5, 17HW E-W MESH S 
MESH SIZE=,F8.5,25HDEG TIME BETWEEN FRAMES^, 


I 38H0TAPE 
I 30H0T I ME 
I 54H0F0LL0WING 
I 13,1 8A4 ) 

I 1 HO , I 8 , 2H , 1 8 A4 


4 FORMAT 

5 FORMAT 

6 FORMAT 

7 FORMAT 

8 FORMAT 
CALL COUNTV 
PRINT 3 

READ (5,7) NTAPE ,T APEN 
PRINT 8 , NTAPE , T APEN 
READ (5,7) NTHR ,THR 
PRINT 8 , NTHR , THR 
READ (5,7) NREPL ,REPL 
PRINT 8,NREPLtREPL 
READ (5,7) NDE V »SDEVN 
PRINT 8, NDEV, SDEVN 


READING AND ARRAY FORMATION TOOK, F10. 5, 7HSEC0NDS 
IN SECONDS FOR OPERATION ,I4,6H WAS ,F15.8) 

ARE TOTALS FOR ALL ARRAY PAIRS IN THIS RUN 


) 


CMXCOOO 1 
CMXC0002 
CMXCOOO- 
CMXC0006 
CMXCOOOS 
CMXCOOOfc 
CMXCOOO 7 
CMXCOOOt 
CMXCOOO 9 
CMXC001C 
CMXC001 1 
CMXC0012 

cmxcooi: 

CMXCOO 1 4 
CMXCOO 13 
CMXCOO 16 
CMXCOO 1 7 
CMXCOO 18 
CMXCOO 1 c . 
CMXC002C 
CMXC002 1 
CMXC0022 
CMXC002 1 
CMXC002^ 
CMXC0025 
CMXC0026 
I CMXC002 7 
CMXC002E 
CMXC002 c . 
CMXC003C 
) CMXC003 1 
CMXC003Z 
) CMXC003 j 
CMXC003 4 
CMXC003 5 
CMXC0036 
CMXC003 7 
CMXC003E 
CMXC003S 
CMXC004C 
CMXC004 1 
CMXC0042 
CMXC0042 
CMXCOO 4^i 
CMXC004 5 
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READ (5,7) NPRN , PRNTN 
PRINT 8 , NPRN , PRNTN 
READ (5,1) MINCLD,MAXCLD, IXZED 
INITIX=IXZED 
I 5 I ZE =88192 
C 

C... 

c 

AMAX=0 
I MAX = 0 
JMAX=0 
I BM=0 
’ JBM=0 
T I M1=0 
TIM2=0 

T I M3=0 ' 

T I M4 = 0 
T I M5=0 

ZI=(0. 0,1.0) 

SET MH ARRAY WITH EXPONENTS OF 2 FOR 64,64,1 TO USE WITH HARM 
\ MH ( 1 ) =6 
> -MH(2)=6 
MH{ 3 ) =0 
MH1=2**MH< 1 ) 

MH2=2**MH( 2 ) 

MHF 1 =MH1 / 2 
MHF2=MH2/2 
MHW 1 =MHF 1+1 
MHW2=MHF2+1 
MWT= MH 1 + 1 
NWT =MH2 + 1 
MPT=MHW1 
NPT=MHW2 
NHARM=0 
1 FS=0 
PRINT 3 

IF (NTAPE.NE.2) GO TO 100 

CALL RDTAPE O 

GO TO 101 

100 READ (5,7) NRN, RN 
PRINT 8 ,NRN , RN 

101 DO 102 I=1,ISIZE 

102 I NT EGR ( I ) =0 

CALL SETHRM ( NHARM , I FE RR , I F S ,MH, MH1 , MH2 , U, I NV, S ) 

CALL TIMEV(RDTIM) 


CMXCOO 46 
CMXC004 7 
CMXC0049 
CMXC0049 
CMXC0050 
CMXC005 1 
CMXC0052 
CMXC005 3 
CMXC005A 
CMXC0055 
CMXCOO 56 
CMXC005 7 
CMXCOO 38 
CMXC005 9 
CMXC0060 
CMXC006 1 
CMXC0062 
CMXC0063 
CMXC0064 
CMXC006 5 
CMXCOO 66 
CMXC0067 
CMXC0068 
CMXC0069 
CMXC0070 
CMXC007 1 
CMXCOO 72 
CMXC007 3 
CMXCOO 74 
CMXC007 5 
CMXC0076 
CMXC007 7 
CMXC0078 
CMXC007 9 
CMXC0080 
CMXC008 1 
CMXCOO 82 
CMXC0083 
CMXC0084 
CMXC008 5 
CMXCOO 86 
CMXC008 7 
CMXC00S8 
CMXC0089 
CMXC0090 
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C PRINT TIME USED TO READ DATA AND SET UP PROGRAM 

PRINT AtRDTIM 
DO 5000 I THRUN= 1 i NTHR 

C INCREMENT COUNTER OF MAXIMUM CLOUD DATA THRESHOLDS USED 

DO 3000 I THRE P= 1 » NRE PL 

C INCREMENT COUNTER OF REPLICATIONS FOR ONE ARRAY PAIR 

READ < 1 5 ) LIST 
IR=LIST(30) 

JR = L I ST ( 31 ) 

) NTOT=IR*JR 

FIRST READ IN THE' NUMBERS OF ARRAY PAIRS E-W AND N-S 
THEN COMPUTE THE TOTAL NUMBER TO BE PROCESSED 


DO 500 I TRY= I tNTOT 
C BEGINNING OF OPERATION 1 
CALL COUNTV 

CALL NULOOP ( U , V , W , T ARA , T ARB , MB, MH 1 , MH2 ) 

READ C 15 ) RL 

READ (15) KARA . 

READ (15) KARB 
X *RINT 2 ,RL 

CALL SETRAY ( RL , DEL I , DEL J , T I MEK ) 

CAL\ ZFRAME < K AR A , MH 1 »MH2 , L AR A ,MHF 1 ,MHF2 , M I NCLD , MAXCLD, NFRQA , 
1PCNTA,, AMEN, SDEVA,VARA»INITIX,, TARA, NDEV) 

CALL ZNORM ( KARB ,MH1 , MH2 , M I NCLD t MAXCLD, NFRQB , PCNTB , BMEN , SDE VB , 
1INITIX»TARB,NDEV) 

CALL T1MNSD ( T ARB , YWMN , YWSD , MH 1 , MH2 ,MHW1 , MHW2 , YMNMN , YMNSD, YSDMN , 
lYSDSDt ZTEST ) 

DO 200 IUV=1 »MH1 
DO 200 JUV=1 1 MH2 

Z( IUV,JUV)=TARA( IUVt JUV)+( ZI*TARB( IUV,JUV) ) 

200 CONTINUE 

C TERMINATION OF OPERATION 1 
CALL TIMEV(TIM( 1, ITRY) ) 

C BEGINNING OF OPERATION 2 
CALL COUNTV 
NOPN= 1 

PRINT 5,N0PN t TIM( NOPN, ITRY) 


1 FS=-2 
NHARM=4 

CALL HARM ( Z t MH , INV , S » I FS , I F ERR ) 

IF ( IFERR.EQ.O) GO TO 201 

CALL IFEXIT ( NHARM ♦ I FERR ♦ I FS t MH, MH 1 ,MH2 » Z » I NV , S ) 
GO TO 5001 


CMXC0091 
CMXC0092 
CMXC009 3 
CMXC0094 
CMXC0095 
CMXC0096 
CMXC009 7 
CMXC0098 
CMXC0099 
CMXC0100 
CMXC0101 
CMXC0102 
CMXC0103 
CMXC0104 
CMXC0105 
CMxcoioe 

CMXC0107 
CMXC0108 
CMXC0109 
CMXC011C 
CMXC0111 
CMXC01 12 

cmxcoii: 

CMXCO 1 14 
CMXC01 1 5 
CMXCO 1 1 6 
CMXC01 1 7 
CMXCO 1 1 E 
CMXC01 1 c . 
CMXCO 1 2( 
CMXC012 3 
CMXCO 122 
CMXCO 12 1 
CMXCO 12^ 
CMXC012! 
CMXC012< 
CMXCO 1 2 ‘ 
CMXCO 121 
CMXC012< 
CMXCO 1 3( 
CMXC013 : 
CMXC013.- 
CMXC013. 
CMXC013- 
CMXC013: 
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201 CONTINUE 

MMNN=MH1*MH2 

CALL SPLITV (MHl,MH2»MMNN,MHFltWtZ J 
TERMINATION OF OPERATION 2 
CALL TIMEV(TIM(2»1TRY) ) 

BEGINNING OF OPERATION 3 
CALL COUNTV . 

N0PN=2 

PRINT 5,N0PN t TIM(N0PN,ITRY) 

I FS=2 
NHARM=3 

CALL HARM '( W»MH, INV»S, IFSt 1FERR) 

IF ( I FERR.EQ.O ) GO TO 211 

CALL IFEXIT ( NHARM t I FERR t I F S t MH , MH 1 t MH2 t W t I NV* S ) 
GO TO 5001 
211 CONTINUE 

TERMINATION OF OPERATION 3 
CALL TIMEV(TIM(3tITRY) ) 

BEGINNING OF OPERATION 4 
^ CALL COUNTV 
N0PN=3 

SPRINT 5 t N0PN,TIM(N0PN» ITRY) 

CALL LAGWS ( W » WT ,MH 1 ,MH2 tMWT ,NWT ) 


GO TO 214 \ 


MxVwCH=0 

IF ( N ZTEST.NE.O) 

MXSWCH= 1 v 

GO TO 215 

214 CONTINUE 

CALL WINDOW ( WT , MWT ,NWT , XCV T MHW 1 ♦ MHW2 » XCVP , MPT , NPT ♦ 
ISDEVAtMXSWCH) 

IF (MXSWCH.NE.O) GO TO 215 

CALL XCMAX ( XCVP , AMAX, MPT , NPT t 1 MAX t JMAX t MXSWCH ) 

215 IF(MXSWCH.EQ.O ) GO TO 220 

CALL WINDOW ( WT , MWT t NWT t XCV ,MHW 1 t MHW2t XCVP » MPT » NPT» 
lSDEVA t MXSWCH) 

IF (MXSWCH.NE.O) GO TO 230 

CALL XCMAX ( XCVP » AM AX ,MPT , NPT , I MAX , JMAX » MXSWCH ) 

IF ( MXSWCH. NE. 6) GO TO 230 

220 CALL X I JMAX ( MPT ,NPT , IM AX t JM AX , I BM t JBM, AMAX » MAXA ) 
CALL VECTOR ( I BM * JBM t DEL I » DEL J, T I MEK , AMAX ) 

C TERMINATION OF OPERATION 4 
230 CONTINUE 

CALL TIMEV(TIM(4,ITRY) ) 

C BEGINNING OF OPERATION 5 
CALL COUNTV 


CMXCO 1 3fc 
CMXCO 1 3 1 
CMXCO 13f 
CMXC013 1 ; 
CMXCO 1 4C 
CMXC01 4 1 
CMXCO 142 
CMXCO 14.' 
CMXCO 1 4* 
CMXCO 14.‘ 
CMXCO 1 4< 
CMXC014' 
CMXCO 14f 
CMXC014 1 
CMXCO 1 5( 
CMXC015. 
CMXC015; 
CMXC015: 
CMXC01 5' 
CMXC015 
CMXC015 
CMXCO 1 5 
CMXCO 15 
CMXC015 
CMXC016! 
CMXC016 
CMXC016 
CMXC016 

YWSDt MHW1,MHW2» CMXCO 16 
CMXC016 
CMXC016 
CMXC016 
CMXC016 

YWSDt MHWl»MHW2t CMXCO 16 
CMXCO 1 7 
CMXC017 
CMXC017 
CMXC017 
CMXCO 1 7 
CMXC017 
CMXCO 1 7 
CMXCO 1 7 
CMXC017 
CMXC017 
CMXCO 1 E 
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ARRAY PRINT-OUTS ARE SPECIFIED 
390 

OF REQUESTED ARRAYS 


N0PN=4 ■ 

PRINT 5,N0PN,TIM(N0PN,ITRY) 

BY-PASS TO 390 IF NO 
IF (NPRN.EQ.O) GO TO 
PROCEED TO PRINT-OUT 
DO 300 K=l,33 
LAGRAY (K,l ) =11111 
300 LAGRAY ( 33 ,K ) =1 1 1 1 1 
DO 305 1=1,32 
DO 305 J=2 ,33 
K = J- 1 ' 

305 LAGRAY ( I »J)=LARA{ I , K ) 

CALL OUTLAGt LAGRAY, 1, IBM, JBM) 

DO 310 K = 1 ,33 
L AGRAY ( K , 1 ) = 1 1 1 1 1 
310 LAGRAY ( 1,K)=11111 
DO 311 K= 1 , 32 
K I = 1 6+K+ I BM 
DO 311 L=2, 33 
LJ=15+L+JBM 

3 Ik LAGRAY ( K , L ) =KARB ( K I , L J ) 

GALL OUTLAG ( LAGRAY ,2 , I BM , JBM ) 

D0\ 320 K= 1 , 33 v 

DO 320 L=1 ,33 \ 

LAGRAY (K,L)=YWMN(K,L) \ 

320 CONTINUE v 

CALL OUTLAGC LAGRAY, 3, IBM, JBM) 

DO 330 K=1 ,33 
DO 330 L= 1 , 33 

LAGRAY ( K,L )=YWSD(K,L)*10.0 
330 CONTINUE 

CALL 0UTLAG(LAGRAY,4, IBM, JBM) 

DO 340 K= 1 , 33 
DO 340 L = 1 ,33 

LAGRAY (K,L)=XCVP(K,L)*400.0 
340 CONTINUE 

CALL OUTLAGI LAGRAY, 5, IBM, JBM) 

TERMINATION OF OPERATION 5 
CALL TIMEV(TIM(5, ITRY) ) 

N0PN=5 

PRINT 5,N0PN,TIM(N0PN, ITRY) 

ARRAY PRINT-OUTS COMPLETED IF SPECIFIED 
390 CONTINUE 
PRINT 3 

ALL PROCESSING FOR THIS ARRAY PAIR IS COMPLETED 


CMXC018 1 

CMXC0182 

CMXC0182 

CMXC018A 

CMXC0185 

CMXC018C 

CMXC0187 

CMXC018E 

CMXCO 18*: 

CMXC019C 

CMXC0191 

CMXC019? 

CMXC0192 

CMXCO 1 9* 

CMXC0195 

CMXC019C 

CMXC0197 

CMXC019? 

CMXC019' 

CMXC020C 

CMXC020 3 

CMXC020; 

CMXC020; 

CMXC020^ 

CMXC020 i 

CMXC020C 

CMXC020' 

CMXC020C 

CMXC020‘ 

CMXC02 1C 

CMXC021 : 

CMXC021; 

CMXC021: 

CMXC021' 

CMXC021 ! 

CMXC02 1< 

CMXC021 

CMXC02 1 1 

CMXC021 1 

CMXC022C 

CMXC022: 

CMXC022. 

CMXC022: 

CMXC022- 

CMXC022! 
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500 

C 


400 

C 

3000 

C 


4000 

C 

5000 


5001 


CONTINUE 

ALL ARRAY PAIRS IN REPLICATION COMPLETED. CUMULATE OPERATIUN 
DO 400 J=1 ,NTOT 
TIM1=TIM1+TIM( 1 , J ) 

TIM2=TIM2+TIM(2,J) 

TIM3=TIM3+TIM( 3, J ) 

TIM 4= T I M 4+T I M ( 4 , J ) * 

TIM5=TIM5+TIM( 5, J ) 

CONTINUE 
REWIND 15 

THIS REPLICATION COMPLETED. RETURN IF MORE SPECIFIED 
CONTINUE 

REPLICATIONS COMPLETED. READ NEW MAX THRESHOLD IF ANY 

NNTH= I THRUN-1 

IF (NNTH.EQ.O) GO TO 4000 

READ (5,7) MAXCLD , AXCLD 

CONTINUE 

REITERATE THROUGH MAIN LOOP IF NEW THRESHOLD READ IN 
CONTINUE 


CMXC0226 
T I ME SCMXC022 7 
CMXC0228 
CMXC022S 
CMXC0230 
CMXC023 1 
CMXC023 2 
CMXC0233 
CMXC0234 
CMXC023 5 
CMXC023E 
CMXC023 7 
CMXC023E 
CMXC023S 
CMXC024C 
CMXC024 1 
CMXC0242 
CMXC0242 
CMXC024A 


PRINT 

3 

CMXC0245 

PRINT 

6 

CMXC0246 

LX= 1 


CMXC0241 

PRINT 

5 »LX ,T I M 1 

CMXC024I 

LX = 2 


CMXC024S 

PRINT 

5»LX,TIM2 

. CMXC025C 

LX = 3 


CMXC025 ] 

PRINT 

5 ,LX,TIM3 

CMXC0252 

LX = 4 


CMXC025; 

PRINT 

5 »LX ,T IM4 

CMXC02 5* 

LX = 5 


CMXC025! 

PRINT 

5,LX,TIM5 

CMXC02 5( 


CONTINUE 

STOP 

END 


CMXC025 ’ 
CMXC0251 
CMXC025* 
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SUBROUTINE BKGRND ( I X , SDEV , AMEAN, V ) 


BKGROOO 1 
BKGR0002 
THAN 10 DBKGROOO 3 
B.KGR0004 
BKGR0005 

AMEAN ARE STANDARD DEVIATION AND MEAN OF PORTION OF SAMPLE BKGR0006 
NOT BACKGROUND (I. E., REEL EC T I V I TY EXCE EDS THRESHHOLD VALUBKGR0007 

BKGR0008 

V IS THE NEXT UNIFORMLY DISTRIBUTED RANDOM VALUE WITH THE SAME MEAN 
STANDARD DEVIATION AS THE NON-BACKGROUND PART OF THIS SAMPLE. 

USE IT TO REPLACE THE BACKGROUND VALUE 


FOR INITIAL ENTRY TO BKGRND, SET IX=AN ODD INTEGER OF LESS 
WE CHOOSE THIS NUMBER TO BE 451798973 


SDE V AND 
WHICH IS 


A = 0. 0 

DO 50 1=1,12 
I Y=IX*65539 
I F ( IY)5,6»6 . 

I Y= I Y+2 1 4748 3647+ 1 

Y = IY 

Y = Y* . 4656613E-9 
I X= I Y 

50 x A=A+Y 

V^(A-6.0)*SDEV+AMEAN 

RETURN 

end\ 


\ 


BKGR0009 
BKGR0010 
BKGR001 1 
BKGROO 1 2 
BKGR001 3 
BKGROO 14 
BKGROO 1 5 
BKGROO 16 
BKGR001 7 
BKGR0018 
BKGR001 9 
BKGR0020 
8KGR002 1 
BKGR0022 
BKGR002 3 
BKGR0024 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 




HARM DISCRETE FOURIER TRANSFORM. BASIC FORTRAN IV 

INPUT PARAMETERS TO BE SET BY USER BEFORE ENTERING HARM- 

A IS A 3-DIMENSIONAL ARRAY OF COMPLEX COE FF I C I ENTS t 
OF DIMENSION! N( I ) ,N( 2 ) ,N( 3) ). 

THE A ' S ARE STORED WITH REAL PART OF A(I1, 12, 13) IN THE LOCATION 
WITH INDEX 2*( I3*N( 1 )*N< 2)+12*N( 1 )+Il ) + l AND THE IMAGINARY PART 
IN THE LOCATION IMMEDIATELY FOLLOWING 

IF THE FOURIER SERIES IS REQUESTED, ARRAY A IS REPLACED BY 
X( J1 , J2 ,J3 )=SUM A!K1,K2,K3)*W1**< K 1* J 1 ) *W2** ( K2* J2 ) * W3** ( K3* J3 ) 
SUMMED OVER K1=0,N<1)-1, K2=0,N(2)-1, K3=0,N(3)-1 

WHERE W I =N ( I ) — TH ROOT OF UNITY. 

MU ) ,1 = 1, 2, 3, WHERE N ( I ) =2**M (I ) IS THE NO. OF PTS.IN THE I-TH.DIM. 
THE DIMENSION OF A IN THE CALLING PROGRAM SHOULD BE TWICE THE 
NUMBER OF COMPLEX ELEMENTS OF THE LARGEST A ARRAY TO BE PROCESS 
ED 

THE COMPLEX X*S ARE STORED IN THE SAME MANNER AS A. 


-.IF THE FOURIER TRANSFORM IS REQUESTED, THE ARGUMENT 


IS TAKEN 


TO BE X 
SERIES. 


AND IS REPLACED BY THE ARRAY A SATISFYING THE FOURIER 


M BEING THE M 


LET MT=MAX(M( 1 ) ,M( 2 ) ,M( 3) )-2, NT=2**MT, WITH 
GIVEN WHEN THE TABLES ARE SET. 

S( J) = SIN( J*PI/!2*NT J ) , J = 1,2, 3, ...NT-1. 

I NV ( J+ 1 ) = WORD CONTAINING BITS OF J IN INVERTED ORDER IN ITS 
RIGHTMOST MT BIT POSITIONS, FOR J = 0 , 1 , 2 , . . . , NT-1 . 

LET I FS=0 TO SET UP SIN AND INV TAB3ES. 

I FS=+1 TO SET UP SIN AND INV TABLES AND DO FOURIER SERIES. 

I FS=-1 TO SET UP SIN AND INV TABLES AND DO FOURIER TRANSFORM 
I FS=+2 TO DO FOURIER SERIES ONLY. 

I FS=-2 TO DO FOURIER TRAN SF0RM o ONLY . 

ONE DOES NOT HAVE TO REPEAT THE CALL TO 'HARM' WITH IFS=0,+1,-1 
IF ONE DOES NOT CHANGE THE MAXIMUM M. 

IFERR=0 IF THE ARGUMENTS M ARE O.K. 

I FERR= 1 IF THERE IS AN ERROR IN CALLING 'HARM' 

IF I FS = 0 » + l , — 1 , IT MEANS THAT THE MAXIMUM M IS GREATER THAN 20 


HARM0001 
HARMOOO 2 
HARM0003 
. HARMOOO A 
HARM0005 
HARMOOO 6 
HARMOOO 7 
HARMOOO 8 
HARM0009 
HARM0010 
HARMOOI 1 
HARM001 2 
HARMOO 1 3 
HARMOOI A 
HARMOO 1 5 
HARMOO 1 6 
-HARMOO 1 7 
HARMOOI 8 
HARMOO 1 9 
HARM002 0 
HARM002 1 
HARM002 2 
HARM0023 
HARM002 A 
HARMOO 2 5 
HARM002 6 
HARM0027 
HARM002 8 
HARM0029 
HARM003 0 
HARM003 1 
HARM003 2 
HARM0033 
HARM003 A 
.HARM0035 
HARM0036 
HARM0037 
HARM003 8 
HARM0039 
HARMOOAO 
HARMOO A 1 
HARMOOA 2 
HARMOO A3 
HARMOOAA 
HARMOO A5 
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10 

12 


13 
1 

14 


16 


18 

20 


OR LESS THAN 3 

IF I FS = +-2 » IT MEANS THAT A SUFFICIENTLY LARGE SIN AND INV TABLE 
HAS NOT BEEN COMPUTED. ONE MUST CALL 'HARM* WITH IFS=0,+-1 AND 
WITH A MAX MCI) GREATER THAN OR EQUAL TO THE MAX MCI) FOR WHICH A 
FOURIER TRANSFORM IS TO BE COMPUTED. 

I FERR = — 1 IF ONE IS CALLING ON ‘HARM*. WITH IFS=0,+-1 TO COMPUTE 
SIN, INV TABLES WHICH IT ALREADY HAS COMPUTED ON A PREVIOUS 
CALL TO HARM WITH THE SAME MAXIMUM M 

REFERENCE- AN ALGORITHM FOR THE MACHINE CALCULATION OF COMPLEX 
FOURIER SERIES, BY J.W. COOLEY AND J.W. TUKEY , MATH. OF COMP. 

VOL. 19, P.297-301, APRIL 1965. 


SUBROUTINE HARM ( A ,M , I NV , S , I F S , IFERR) 

DIMENSION AC 1 ) , INVC 1 ) ,S< 1 ) ,N( 3) ,M< 3) ,NP (3 ),W<2), W2C2), W3C 2) 
EQUIVALENCE ( N 1 ,N ( 1 ) ) , ( N2 ,N ( 2 ) ) , ( N3 ,N C 3 ) ) 

IF ( IABSC IFS )-l ) 900,900,12 
MTT=MAXO(M( 1 ) ,M(2) ,M{ 3) ) -2 
ROOT 2 = SORT ( 2 . ) 

IF CMTT-MT ) 14,14,13 
I FERR= 1 
RETURN 
I FERR = 0 
Ml =M ( 1 ) 

M2 = M ( 2 ) 

M3 = M ( 3 ) 

N1=2**M1 
N2=2**M2 
N3=2**M3 
IF (IFS) 16,1,20 

TO CALCULATE TRANSFORM REPLACE A BY CONJG(A)/N 

NTOT = N1*N2*N3 

FN = NTOT 

DO 18 1=1, NTOT 

AC 2*1 — 1 ) = AC 2*1-1 ) / FN 

A ( 2 * I. ) = -AC 2*1 ) /FN 

NPC 1 ) =N1 *2 

NP ( 2 ) = NPC 1 )*N2 

NP ( 3 ) =NP ( 2 ) *N3 

DO 250 10=1,3 

IL = NP ( 3 )-NP ( ID) 

I L 1 = IL + 1 
MI = M(.ID) 


HARM0046 
HARM0047 
HARM0048 
HARM0049 
HARM0050 
HARM005 1 
HARM005 2 
HARM0053 
HARM005 4 
HARM0055 
HARM005 6 
HARM0057 
HARM0058 
HARM0059 
HARM0060 
HARM0061 
HARM0062 
HARM0063 
HARM0064 
HARM0065 
HARM0066 
HARM0067 
HARM0068 
HARM0069 
HARM0070 
HARM007 1 
HARM0072 
HARM0073 
HARM0074 
HARM0075 
HARM007 6 
HARM0077 
HARM0078 
HARMOO 79 
HARM0080 
HARM0081 
HARM0082 
HARM0083 
HARM008 4 
HARM0085 
HARM008 6 
HARM0087 
HARM0088 
HARM0089 
HARM0090 
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IF (MI ) 250 » 2 50 » 30 
30 I D I F = NP ( ID) 

KB 1 T = NP ( ID) 

MEV = 2*(MI/2) 

IF (MI - MEV )60,60,40 
C M IS ODD. DO L=1 CASE 

40 ■ KB I T = KB IT/2 
KL=KB I T-2 

DO 50 I=1,IL1,IDIF 
KLAST=KL+ I 
DO 50 K = I tKLAST ,2 
KD=K+KB I T 

DO ONE STEP WITH L=1,J=0 
A ( K ) = A ( K ) +A ( KD ) 

A ( KD ) =A ( K ) -A { KD ) 

T=A ( KD ) 

A{ KD )=A( K )— T 
A ( K ) =A ( K ) +T 
T=A( KD+1 ) 

A ( KD+1 ) =A ( K+l ) -T 
50 A ( K+ 1 ) =A ( K+l ) +T 

IF (MI - 1)250,250,52 
52 LFIRST =3 

DEF - J LAST = 2**(L-2) -1 
JLAST= 1 
GO TO 70 
M IS EVEN 
60 LFIRST = 2 
JLAST=0 

70 DO 240 L=LFIRST,MI ,2 
JJDIF=KBIT 
KBIT =KB I T/4 
KL=KB I T-2 
DO FOR J=0 
DO 80 1 = 1,1 LI, IDIF 
KLAST= I +KL 
DO 80 K= I , KLAST , 2 
K1=K+KB IT 
K2=K 1+KB I T 
K3=K2+KB I T 

DO TWO STEPS WITH J=0 
A(K)=A(K)+A(K2) 

A ( K2 ) = A ( K ) —A ( K2 ) 


HARM0091 
HARM0092 
HARM0093 
HARM0094 
HARM0095 
HARM0096 
HARM0097 
HARM0098 
HARM0099 
HARM0100 
HARM0101 
HARM0102 
HARM0103 
HARM0104 
HARM0105 
HARM0106 
HARM0107 
HARM0108 
HARM0109 
HARM0110 
HARMO 111 
HARM01 1 2 
HARMO 113 
HARM0114 
HARMO 115 
HARM01 1 6 
HARMO 117 
HARM01 1 8 
HARM01 19 
HARM0120 
HARM0121 
HARM0122 

HARM0123 

HARM0124 
HARM0125 
HARM0126 
HARMO 127 
HARM012 8 
HARMO 129 
HARM0130 
HARM0131 
HARM01 32 
HARMO 133 
HARM0134 
HARM01 35 
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A ( K 1 ) = A ( K 1 ) +A ( K 3 ) 
A(K3)=A(K1 )— A(K3 ) 

A ( K ) =A ( K ) +A ( K 1 ) 

A ( K 1 ) = A ( K )- A ( K 1 ) 

A(K2 ) = A ( K2 )+A(K3 )*I 
A ( K3 ) = A ( K2 ) -A ( K3 ) * I 

T=A(K2) 

A ( K2 ) =A ( K ) — T 
A ( K ) =A ( K ) +T 
. T=A ( K2+1 ) 

A ( K2+ 1 ) =A ( K+l ) -T 
A(K+1 )=A(K+1 )+T 

T=A ( K3 ) 

A ( K3 ) = A ( K 1 )-T 
A ( K 1 ) =A ( K 1 ) +T 
T= A ( K3+ 1 ) 

■-) A ( K 3+ 1 )=A(K1 + 1 )— T 
\_.A(K1 + 1 )=A(K1+1 )+T 

T = A( K1 ) 

A ( K I ) =A ( K ) -T 
A(K)=A(K )+T 
T=A ( K 1 + 1 ) 

A ( K 1 + 1 )=A(K + 1 )-T 
A(K+1 ) =A ( K+l )+T 

R=— A ( K3+1 ) 

T = A ( K3 ) 

A(K3)=A(K2)-R 
A ( K2 ) = A ( K2 ) +R 
A( K3+1 ) = A( K2 + 1 )— T 
80 A ( K2+ 1 ) =A { K2+ 1 ) +T 

IF ( JLAST ) 235,235,82 
82 J J= J JD I F +1 

DO FOR J=1 

I LAST= IL +JJ 

DO 85 I = JJ, HAST, IDIF 

KLAST = KL+I 

DO 85 K= I , KLAST , 2 

K1 = K+KBIT 

K2 = KI+KBIT 


HARM0136 
HARM0I37 
HARMOI38 
HARM0135 
HARMOIAC 
HARM01A1 
HARM0142 
HARM01A3 
HARM01AA 
HARM0145 
HARM01 A 6 
HARM0147 
HARMOlAc 
HARM014S 

HARM015C 
HARM0151 
HARMO 152 
HARM0153 
HARM015^ 
HARMO 155 
HARMO 156 
HARMO 15*. 
HARM015J 
HARMO 1 5 S 
HARM016C 
HARMO 1 63 
HARM016; 
HARM0162 
HARM016' 
HARMO 1 6f 
HARMO 1 6< 
HARMO 1 6' 
HARM016f 
HARMO 1 6‘ 
HARMO 1 7 ( 
HARM017: 
HARMO 17. 
HARMO 1 73 
HARM017* 
HARMO 1 7: 
HARMO 17 ; 
HARM017 
HARM017 1 
HARMO 1 7‘ 
HARMO 1 8 < 
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K3 = K2+KBIT 

LETTING W=< 1 + 1 )/ROOT2,W3=(-l+I ) /R00T2, W2=I , 

A(K )=A( K )+A(K2 )*I 

A<K2)=A(K)-A(K2)*I 

A(K1)=A(K1)*W+A<K3)*W3 

A(K3)=A(K1 )*W-A(K3)*W3 

A ( K ) = A ( K ) +A ( K 1 ) 

A ( K 1 ) =A ( K ) —A ( K 1 ) 

A(K2)=A(K2)+A(K3)*I 
A(K3)=A(K2 ) —A ( K3 ) * I 

R =— A ( K2 + 1 ) 

T = A ( K2 ) 

A ( K2 ) = AIK) — R 
AIK) = AIKJ+R 
A ( K2+ 1 ) =A ( K+ 1 ) — T 
A I K+ 1 ) =A ( K + l ) +T 

AWR=A(K1 ) — A (Kl + 1 ) 

AW I = A (Kl + 1 )+A{ K1 ) 

R=-AIK3)— A(K3+1 ) 

T=A ( K3 ) —A I K3 + 1 ) 

A ( K3 ) = ( A WR-R ) /ROOT 2 
AIK3+1 )=( AWI-T 5/ROOT2 
A(K1 )=(AWR+R)/R00T2 
AIKl+l )=(AWI+T )/R00T2 
T = A ( K 1 ) 

A ( K1 ) =A ( K )— T 
A ( K ) =A ( K ) +T 
T=A ( K 1+1 ) 

A ( Kl+1 )=A ( K+l )— T 
A ( K + l ) =A ( K + l ) +T 
R=— A ( K3+1 ) 

T=A(K3) 

A ( K3 ) = A ( K2 )-R 
A ( K2 ) =A ( K2 ) +R 
AIK3+1 ) = A ( K2+ 1 )-T 
85 AIK2+1 )=A(K2+1 )+T 

I F ( JLAST-1 ) 235 » 235 » 90 
90 JJ= JJ + JJDIF 

NOW DO THE REMAINING J'S 
DO 230 J=2»JLAST 


HARM0181 
HARM0182 
HARM0183 
HARM018A 
HARM0185 
HARM0186 
HARM0187 
HARM0188 
HARM0189 
HARM0190 
HARM0191 
HARM0192 
HARM01 93 
HARM0194 
HARM0195 
HARM0196 
HARM0197 
HARM0198 
HARM0199 
HARM0200 
HARMQ.2.01 
. HARM0202 
HARM0203 
HARM0204 
HARM0205 
HARM0206 
HARM0207 
HARM0208 
HARM0209 
HARM021 0 
HARM02 1 1 
HARM021 2 
HARM02 13 
HARM021 4 
HARM0215 
HARM021 6 
HARM02 17 
HARM021 8 
HARM02 19 
HARM0220 
HARM022 1 
HARM0222 
HARM0223 
HARM0224 
HARM0225 
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b f TCH W* S 

Off:- w=W**INV( J ) , W2=W**2, W3=W**3 
J = INV( J+l ) 

<.<- 1 ( = NT - 1 
»< 1 >.= S< IC) 
w I 2 > = S ( I ) 

12 = 2*1- 
12C=NT-12 

JF< 120120,110,100 

2* I IS IN FIRST QUADRANT 
t on w2( 1 )=S( I 20 
w 2 ( 2 ) = S ( 12) 

GO TO 130 

no w 2 ( l )=o. 

w 2 ( 2 ) = 1 . 

GO TO 130 


(, 2*1 IS IN SECOND QUADRANT 

-4 2 0 1 2CC = I2C+NT 
\ 12C=-I2C 

X?<1 ) =-S ( I 2 C ) 
w>12) = S( I2CC) 
no I 3 =\+ 1 2 
I 3C = NT- 1 3 

1H 130160, 150,140 


C 13 IN FIRST QUADRANT 

1^0 W3t 1 ) =S C I3C) 

W3(2)=S( 13) 

GO TO 200 
l*>o w3( 1 } =0 . 

H3(2)=l. 

GO TO 200 


no I 3CC= 1 3C + NT 

IF( I 3CC ) 190 t 180, 170 


I 3 IN SECOND QUADRANT 
1 tO I 3C=- I 3C 

W3( 1 J=-S( I3C J 
*3 ( 2 ) =S ( I3CC) 

GO TO 200 
no w3i l n-i. 

*<3(2 ) = 0. 


HARM022 6 
HARM0227 
HARM022 8 
HARM0229 
HARM0230 
HARM0231 
HARM0232 
HARM0233 
HARM0234 
HARM0235 
HARM0236 
HARM0237 
HARM0238 
HARM0239 
HARM0240 
HARM0241 
HARM0242 
HARM0243 
HARM0244 
HARM0245 
HARM0246 
HARM0247 
HARM0248 
HARM0249 
HARM0250 
HARM02 5 1 
HARM0252 
HARM0253 
HARM0254 
HARM02 55 
HARM025 6 
HARM02 57 
HARM0258 
HARM02 59 
HARM0260 
HARM0261 
HARM0262 
HARM0263 
HARM0264 
HARM0265 
HARM0266 
HARM0267 
HARM0268 
HARM0269 
HARM0270 
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GO TO 200 
C 

C 3*1 IN THIRD QUADRANT 

190 I 3CCC = NT + 1 3CC 
I 3CC • = -I3CC 
W3( 1 )=-S( I3CCC) 

W3 I 2 ) =-S I I 3CC ) 

200 I L AST= I L + J J 

DO 220 I=jj»ILAST»IDIF 

KL AST=KL+ I 

DO 220 K = I tKLAST »2 

K 1 =K+KB I T 

K2=K 1 +KB I T 

K3=K2+KB I T 


DO TWO STEPS WITH J NOT 0 
AIK )=A IK )+A(K2 )*W2 
AIK2)=AIK)-A(K2)*W2 
A ( K 1 ) =A ( K 1 )*W+A(K3)*W3 
A I K3 ) = A I K 1 ) * W— A ( K 3 ) * W3 


2WK)=A(K)+A(K1) 

A (Xl )=A(K )-A(Kl ) 
A(K2 n )=A(K2)+A(K3)*1 
A(K3)=A(K2 J-AIK3 )*I 



R = A(K2)*W2( 1 ) — A ( K2 + 1 )*W2( 2 ) 

T=A(K2)*W2(2)+A(K2+1 )*W2( 1) 

A ( K2 ) =A ( K ) — R 

AIK) =A ( K ) +R 

AIK2+1 )=A(K+1 ) -T 

AIK+l )=A( K + l-J+T 

R = A(K3)*W3( 1 ) -A I K3 + 1 )*W3I2) 
T=AIK3)*W3(2)+AIK3+1)*W3( 1) 
AWR=A I K1 )*W( 1 J-AIKl+l )*W< 2) 
AWI=A(K1 )*W(2 )+A(Kl+l)*WI 1) 
A ( K3 ) = AWR— R 
A I K3+1 > =AW I — T 
A I K 1 ) = AWR+R 
AIKl+l )=AWI+T 
T = A I K 1 ) 

A I K 1 ) =A I K ) — T 
A I K ) =A I K ) +T 
T=A I K 1+ 1 ) 


HARM0271 

HARM0272 

HARM0273 

HARM0274 

HARM02 75 

HARM0276 

HARM02 77 

HARM0278 

HARM02 79 

HARM0280 

HARM0281 

HARM028 2 

HARM0283 

HARM0284 

HARM0285 

HARM0286 

HARM02 87 

HARM0288 

HARM0289 

HARM0290 

HARM0291 

HARM0292 

HARM0293 

HARM0294 

HARM0295 

HARM0296 

HARH0297 

HARM029E 

HARM029S 

HARM030C 

HARM0301 

HARM0302 

HARM0302 

HARM030* 

HARM0305 

HARM030d 

HARM030' 

HARM030I 

HARM030C 

HARM031C 

HARM03 1 I 

HARM031 •' 

HARM031; 

HARM031 / 

HARM03 1 t 
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A(K1+1 )=A(K+1 )— T 
A ( K+l ) = A ( K+ 1 ) +T 
R=-A ( K3+ 1 ) 

T=A(K3) 

A( K3)=A( K2 )-R 
A ( K2 ) =A ( K2 ) +R 
A(K3+1 )=A(K2+1 )-T 
220 A<K2+1 )=A(K2+1 ) +T 
C END OF I AND K LOOPS 

230 J J= J JDI F+J J 
C END OF J-LOOP 

235 JLAST=4*JLAST+3 
240 CONTINUE 
C END OF L LOOP 

250 CONTINUE 

END OF ID LOOP 

WE NOW HAVE THE COMPLEX FOURIER SUMS BUT THEIR ADDRESSES ARE 
BIT-REVERSED. THE FOLLOWING ROUTINE PUTS THEM IN ORDER 
N NTSQ=NT*NT 
-,M3MT = M3-MT 

350 I F ( M3MT ) 370 t 360 T 360 
C M3 GR. OR EQ. MT 

360 I G03= 1 

N3VNT=N3/NT 
M I NN3=NT 
GO TO 380 

C M3 LESS THAN MT 

370 I G03=2 
N3VNT = 1 
NTVN3=NT/N3 
MINN3=N3 

380 JJD3 = NTSQ/N3 
M2MT=M2-MT 

450 IF { M2MT )470 t 460,460 
C M2 GR. OR EQ. MT 

460 IG02=1 ^ 

N2VNT =N2/NT 
M I NN2=NT 
GO TO 480 

C M2 LESS THAN MT 

470 IG02 =2 
N2VNT = 1 
NTVN2 =NT/N2 
M I NN2=N2 


HARM03 1 6 
HARM0317 
HARM03 1 8 
HARM03 1 9 
HARM032 0 
HARM03 2 1 
HARM032 2 
HARM0323 
HARM0324 
HARM0325 
HARM032 6 
HARM0327 
HARM032 8 
HARM0329 
HARM0330 
HARM0331 
HARM0332 
HARM033 3 
HARM0334 
HARM0335 
HARM0336 
HARM0337 
HARM03 3 8 
HARM0339 
HARM0340 
HARM0341 
HARM0342 
HARM0343 
HARM0344 
HARM0345 
HARM0346 
HARM0347 
HARM0348 
HARM0349 
HARM035 0 
HARM0351 
HARM035 2 
HARM03 53 
HARM0354 
HARM03 55 
HARM035 6 
HARM0357 
HARM0358 
HARM03 59 
HARM0360 


F-17 



) 


480 J JD2 = Nt SQ/N2 
M1MT=M1-MT 

550 IF(M1MT >570,560,560 
C Ml GR« OR EO. MT 

560 I GO 1 = 1 

N1VNT=N1/NT 
M I NN 1 =NT 
GO TO 580 

C Ml LESS THAN MT 

570 IG01=2 
N 1 VNT = 1 
NTVN 1=NT /N1 
MINN1=N1 

580 JJD1=NTS0/N1 
600 J J3=l 
J=1 


DO 880 JPP3 = 1 , N3VNT 
IPP3=INV( JJ3) 

DO 870 JP3=1 »M INN3 
\ GO TO (610,620), I G03 
610\IP3=1NV ( JP3 ) *N3VNT 
G^ TO 630 

620 1P3=INV( JP3 )/NTVN3 
630 I3=f\JPP3+IP3)*N2 
700 JJ2=1 


00 870 JPP2=1,N2VNT 
1PP2=INV( JJ2 )+I3 
DO 860 JP2=1 ,MINN2 
GO TO (710, 720), I G02 
710 IP2=INV< JP2)*N2VNT 
GO TO 730 

720 IP2=INV(JP2)/NTVN2 
730 1 2= ( I PP2 + 1 P2 ) *N 1 
800 JJ1=1 


DO 860 JPP1 =1 »N1VNT 
IPP1=INV( JJ1 ) + 1 2 
DO 850 JP1 = 1, MINN-1 
GO TO (810,820) , IG01 
810 I P 1 = I NV ( J P 1 )*N1VNT 
GO TO 830 

820 IP1=INV( JP1 J/NTVN1 
830 I=2*( IPP1+IP1 )+l 

IF ( J-I ) 840,845,845 
840 T = A ( I ) 

A ( I ) = A ( J ) 



i. 



HARM0361 

HARM036 2 

HARM0363 

FIARM0364 

FiARM0365 

HARM0366 

HARM0367 

HARM0368 

HARM0369 

HARM0370 

HARM0371 

HARM037 2 

HARM0373 

HARM037 4 

HARM03 75 

HARM0376 

HARM0377 

HARM0378 

HARM03 79 

HARM038C 

HARM0381 

HARM0382 

HARM0383 

HARM0384 

HARM0385 

HARM038 6 

HARM03 87 

HARM038E 

HARM038S 

HARM039C 

HARM0391 

HARM0397 

HARM0393 

HARM039* 

HARM0395 

HARM039( 

HARM039', 

HARM039F 

HARM0399 

HARM040C 

HARM040 ] 

HARM040; 

HARM040; 

HARM040' 

HARM0405 
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A(J)=T 
T=A ( 1 +1 ) 

A( 1 + 1 )=A( J + l ) 

A ( J + l ) =T 
845 CONTINUE 
850 J= J+2 
860 JJ1=JJ1+JJD1 
C END OF JPPI AND JP2 

870 JJ2= JJ2+JJD2 
C END OF JPP2 AND JP3 LOOPS 

880 JJ3 = JJ3+JJD3 ' 

C END OF JPP3 LOOP 

' IF(IFS) 882*1*1 

DOING TRANSFORM. REPLACE A BY CONJG(A). 

882 DO 884 I=1,NT0T 
884 A { 2* I ) = -A (2*1 ) 

GO TO 1 
RETURN 

THE FOLLOWING PROGRAM COMPUTES THE SIN AND INV TABLES. 

900 MT=MAX0(M(1)*M(2)*M(3) ) -2 
MT = MAXO ( 2 »MT ) 

904 IF <MT-20>906,906,905 . 

905 IFERR =1 
GO TO 1 
RETURN 

906 I FERR=0 ’ 

NT=2**MT 
NTV2=NT/2 
SET UP SIN TABLE 
THETA=PIE/2**(L+1 ) FOR L=1 

910 THETA=. 7853981634 
C JSTEP=2**(MT-L+1 ) FOR L=1 

JSTEP=NT 

C JDIF=2**(MT-L) FOR L=1 

JDI F=NT V2 ‘3 

S(JDIF)=SIN(THETA) 

DO 950 L=2 *MT 
THETA=THETA/2. 

JSTEP2= JSTEP 
JSTE P = JD I F 
JD I F- JSTEP/2 
S(JDIF)=SIN(THETA) 


HARM040 6 
HARM0407 
HARM040 8 
HARM0409 
HARK0410 
HARM04 1 1 
HARM041 2 
HARM04 1 3 
HARM0414 
HARM04 1 5 
HARM04 1 6 
HARM04 1 7 
HARM04 1 8 
HARM0419 
HARM042 0 
HARM042 1 
HARM042 2 
HARM0423 
HARM042 4 
HARM0425 
HARM042 6 
HARM0427 
HARM042 8 
HARM0429 
HARM043 0 
HARM043 1 
HARM0432 
HARM0433 
HARM0434 
HARM0435 
HARM043 6 
HARM0437 
HARM0438 
HARM0439 
HARM0440 
HARM0441 
HARM0442 
HARM0443 
HARM0444 
HARM0445 
HARM0446 
HARM0447 
HARM0448 
HARM0449 
HARM0450 
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JC 1 =NT- JD I F 
S< JC1 )=COS< THETA ) 
JLAST=NT-JSTEP2 
I F ( JLAST - JSTEP) 950,920,920 
920 DO 940 J=JSTEP, JLAST, JSTEP 
JC=NT- J 
JD = J + JD IF > 

940 S(JD)=S(J)*S(JC1 )+S(JDIF)*S(JC) 
950 CONTINUE 

SET UP INV(J) TABLE 

960 MTL EXP = NT V2 
C MTLEXP=2**(MT-L ). FOR L=1 

LM1 EXP = 1 

C LM1EXP=2**(L-1). FOR L=1 

I NV ( 1 ) =0 
DO 980 L = 1 , MT 
INVILM1EXP + 1 ) = MTLEXP 
DO 970 J=2,LM1EXP 
J J = J + LM 1 EXP 

970 INV( JJ ) - I NV { J I+MTLEXP 
MTLEXP=MTLEXP/2 
980 LM1EXP=LM1EXP*2 
I F ( IFS) 12,1,12 
C RETURN 

END 


HARM04 5 1 
HARM045 2 
HARM0452 
HARM0454 
HARM0455 
HARM045 6 
HARM0457 
HARM045E 
HARM0459 
HARM046C 
HARM0461 
HARM0462 
HARM0462 
HARM046< 
HARM0465 
HARM0466 
HARM0467 
HARM046E 
HARM0469 
HARM047C 
HARM0471 
HARM0472 
HARM0472 
HARM0474 
HARM0475 
HARM047 6 
HARM0477 
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SUBROUTINE IFEXIT ( NHARM r I FERR , I F S » M , MM, NN , U, I NV , S ) I 

DIMENSION M( 3 ) ,U(MM,NN ) , INV(MM) , S( MM ) I 

COMPLEX U I 

6 FORMAT ( 4 9 HO AN ERROR WAS DETECTED BY HARM AT ENTRY NUMBER ,13, I 

112H IFERR WAS ,I3,10H IFS WAS ,13 ) I 

16 FORMAT ( 1 1 9H0PR0GRAM HAS BEEN RESET FOR A 64X64 COMPLEX ARRAY AND I 
1A SECOND ATTEMPT IS BEING MADE TO I NI T I AL I Z E SIN AND INV TABLES. )I 
26 FORMAT ( 54H0 SECOND TRY WITH 64X64 ARRAY ALSO FAILED. JOB STOPPED) I 

PRINT 6, NHARM, IFERR, IFS I 

1 FERR=0 I 

IF (NHARM. GE.l) GO TO 1 I 

C NHARM=0 SO ERROR OCCURRED IN TRYING TO INITIALIZE INV AND SIN TABLES I 
C PROGRAM WILL SET HARM FOR 64X64 COMPLEX ARRAY AND PROCEED I 

PRINT 16 I 

M ( 1 ) =6 1 

M( 2 ) =6 'I 

M ( 3 ) =0 I 

CALL HARM ( U, M , I NV , S , 0 , I FERR ) * I 

IF ( IFERR. EO.O) GO TO 3 • I 

' "N PRINT 26 1 

H. CONTINUE I 

NHARM =99 1 

3 CONTINUE I 

RETURN I 

END I 


FEX00U1 
FEX0002 
FEX0003 
FEX0004 
FEX0005 
FEX0006 
FEX0007 
FEX0008 
FEX0009 
FEX0010 
FEX001 1 
FEX0012 
FEX0013 
FEX0014 
FEX0015 
FEX0016 
FEX0017 
FEX0018 
FEX0019 
FEX0020 
FEX0021 
FEX0022 
FEX0023 
FEX0024 
FEX0025 


* 0 . 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE LAGWS I WS , WT , MWS , NW S , MWT , NWT ) LAGWOOOl 

lagwooo; 

THIS REORGANIZES DATA OF COMPLEX ARRAY WS(MWStNWS) AND GENERATES A LAGWOOO; 
NEW ARRAY WT(MWT,NWT) WHERE MWT=MWS+1 AND NWT=NWS+1 WHICH CONTAINS THLAGWOOO* 
SAME DATA AS IT WOULD APPEAR IN LAG COORDINATES. FOR MWS=64 AND NWS=6LAGW000^ 


THE RANGE OF I LAGS AND JL AGS IN WT WIL1 BE FROM -32 TO +32. THE POS I TLAGW0006 


OF DATA CORRESPONDING TO 

I L AG= - 

32, -16, 

-08, 

0 , +0 8, +16, +32 

LAGWOOOl 

WIL HAVE THE COORDINATE 

1 = 

01, 17, 

25, 

33, 41, 49, 65 

LAGWOOOf 

AND 

OF DATA CORRESPONDING TO 

JL AG = 

+32, +16, 

+08, 

0, -08, -16, -32 

LAGWOOOL 

LAGW001C 

WILL HAVE THE COORDINATE 

J = 

65, 49, 

41, 

33, 25, 17, 01 

LAGW001] 

DIMENSION WS ( M WS , NWS ) , 
COMPLEX WS , WT 
KIT=MWS/2 

WT ( MWT 

, NWT ) 



lagwooi; 

lagwooi: 

LAGWOO 1 l 
LAGWOO 1 l 

LIT =K I T+ 1 
M I T = MWS + 1 





LAGWOOK 

LAGWOOI' 

K I S=K I T 
L I S=K 1 S 
K JT=NWS/2+l 




- 

LAGWOOlf 

LAGWOOI* 

LAGW002( 


L JT = K JT + 1 LAGWOO'2: 


MJT=NWS+ 1 
K JS=K JT+1 
L JS=NWS+KJS 
DO 30 I T=1 »K IT 
I S= I T+KI S 
DO 10 JT=1,KJT 
JS=KJS-JT 

WT { IT , JT )=WS( ISrJS ) 
10 CONTINUE 

DO 20 JT=LJT,MJT 
JS=L JS-JT 

WT ( I T , JT ) =WS ( I S , JS ) 
20 CONTINUE 
30 CONTINUE 

DO 60 IT=LIT,MIT 
I S= IT-L I S 
DO 40 JT=1,KJT 
JS=KJS— JT 

WT( IT, JT )=WS( IS, JS) 
40 CONTINUE 

DO 50 JT = LJT ,MJT 
JS=L J S— JT 

WT ( I T , JT ) = WS(ISyJS) 
50 CONTINUE 


LAGW002i 

LAGW002: 

LAGW002 / 

LAGW0021 

LAGW002 < 

LAGW002' 

LAGW002 ! 

LAGW002 1 

LAGW003( 

LAGW003 : 

LAGW003. 

LAGW003: 

LAGW003 

LAGW003! 

LAGW003 

LAGW003 

LAGW003! 

LAGW003' 

LAGW004( 

LAGW004 

LAGW004 

LAGW004 

LAGW004 

LAGW004 1 
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60 CONTINUE 
RETURN 
END 


LAGWOOAf 

LAGW0041 

LAGW004* 



SUBROUTINE NULOOP ( U » V , W, A » B T MS t MM t NN ) 

DIMENSION U( MM *NN ) ♦ V( MM ,NN ) f W( MM»NN) » A( MM»NN ) t B< MM»NN ) »MS( MM»NN ) 
COMPLEX U»V,W 
DO 20 I =1 , MM 
DO 20 J= 1 » NN 
U ( I * J ) =0 

V ( I t J ) =0 . 

W( I » J )=0 
A ( I » J ) =0 


B ( I » J ) =0 
MS ( I » J ) =0 
CONTINUE 
RETURN 
END 


NULOOOO 1 
NUL00002 
NULOOOO 3 
NUL00004 
NULOOOO 5 
NUL00006 
NUL00007 
NUL00008 
NUL00009 
NUL00C10 
NUL0001 1 
NULOOO 12 
NULOOO 1 3 
NULOOO 1 A 



ooooooooooooonnoooonooooo 


SUBROUTINE OUTLAG < LAGRAY ,NDATA, IBM, JBM) 0UTL0001 

OUT L000 2 
0UTL0003 

THIS PRINTS OUT THE CONTENTS OF AN ARRAY ORGANIZED IN LAG FORMAT WITH0UTL0004 
LAGS RANGING FROM -16 TO +16 IN THE ARRAY AND WITH THE HEADER FORMAT 0UTL0005 
NUMBERED NDATA. DATA IN THE ARRAY MUST BE IN INTEGER FORMAT AND MUST 0UTL0006 
NUT EXCEED 999 IN VALUE. ONLY THOSE ELEMENTS CORRESPONDING TO I LAGS 0UTL0007 
C'F -16 THROUGH +15 WILL BE PRINTED ACROSS THE PAGE , THOUGH VERT I CAL L YOU TL OOO 8 
ALL ROWS CORRESPONDING TO I L AGS OF +16 DOWN THROUGH 


-16 WILL APPEAR 


IF ONE ELEMENT OF THE ARRAY HAS SPECIAL MEANING 
MUST BE INCLUDED IN THE CALL AS IBM, JBM. 


ITS LAG COORDINATES 


FOR NDATA = 1 PRINTS DATA VALUES OF WINDOW ARRAY (AXES TRANSLATED) 


PRINTS 

BEST 


DATA VALUES 
MATCH (WITH 


OF SEARCH AREA SUBARRAY SELECTED 
AXES TRANSLATED) 


3 PRINTS MEAN VALUES OF SUBARRAYS OF SEARCH AREA BY 

4 PRINTS STANDARD DEVIATIONS OF SUBARRAYS OF SEARCH 


5 PRINTS CROSS-CORRELATIONS IN PERCENTAGES FOR ALL LAGS 


DIMENSION LAGRAY(33»33) »ILAG(32),JLAG(33) 

1 FORMAT ( 86H1DATA USED AS WINDOW ARRAY WITH COORDINATE AXES MOVED 
ISO COORDINATES APPEAR AS LAGS ) 

2 FORMAT ( 1 2 BH ID AT A CONTAINED IN THE SEARCH 

1 AS THE BEST MATCH. COORDINATE AXES MOVED 

2 LAGS ) 

3 FORMAT ( 1 16H1TRUNCATED MEANS OF SUBARRAYS OF SEARCH AREA 
1DING TO THE LAG VALUES OF WINDOW AREA RELATIVE TO SEARCH 


AREA SUBARRAY 
SO COORDINATES 


OUTL0009 
0UTL0010 
OUT LOO 1 1 
OUTLOO 1 2 
OUTLOO 13 
0UTL0014 
OUTLOO 1 5 
0UTL001 6 
AS0UTL0017 
0UTL001 8 
0UTL0019 
LAGS0UTL0020 
OUT L002 1 
AREA0UTL0022 
0UTL0023 
OUTL0024 
0UTL0025 
0UTL002 6 
0UTL0027 
0UTL0028 
0UTL0029 
OUTL0030 
IDENTIFIED0UTL0031 
APPEAR AS0UTL0032 
0UTL0033 
CORRESPONOUTL0034 
AREA ) 0UTL0035 


4 FORMAT ( 1 1 9H1 STANDARD DEVIATION (TENTHS) OF SEARCH AREA SUBARRAY COUTL0036 
20RRESP0NDI NG TO LAG VALUES OF WINDOW AREA RELATIVE TO SEARCH AREA )0UTL0037 

5 FORMAT ( 129H1CR0SS-C0RRELATI0NS (IN PERCENTAGE) BETWEEN WINDOW AND0UTL0038 
1 SEARCH AREA SUBARRAY ARRANGED BY LAG COORDINATES RELATIVE TO SEAR0UTL0039 


2CH AREA ) 

10 FORMAT (5H0LAGS,I3,31I4> 

11 FORMAT ( 1H0, 13,3214) 

DO 100 1=1,32 

I LAG ( I ) = I - 1 7 
JLAG( I ) = I -1 7 


OUTL0040 

0UTL0041 

0UTL0042 

0UTL0043 

0UTL0044 

0UTL0045 
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100 CONTINUE 

JLAGt 33 ) =1 6 

GO TO (21,22,23,24,25),NDATA 
GO TO 30 

21 PRINT 1 
GO TO 30 

22 PRINT 2 
GO TO 30 

23 PRINT 3 
GO TO 30 

24 PRINT 4 
GO TO 30 

25 PRINT 5 
30 CONTINUE 

PRINT 10, ( ILAG( I ), 1 = 1,32) 

DO 200 1=1,33 
J=34-I 

PRINT 11, JLAG(J)»(LAGRAY(K,J)»K=1» 
200 CONTINUE 
RETURN 
END 


0UTL0046 
OUTL0047 
OUT L004E 
OUT L0049 
0UTL005C 
OUTL0051 
0UTL0052 
0UTL0053 
0UTL0054 
0UTL0055 
0UTL005C 
0UTL0057 
0UTL005E 
OUT LOO 59 
0UTL006C 
0UTL0061 
OUTL0062 

) < OUTL0063 

- OUTL006^ 

0UTL0065 
OUTL006E 



SUBROUTINE RAYSET ( L t R t 1 1 J « P t Nt IR't JR ) 
DIMENSION R ( I R T JR 1 5 ) ♦ L < 32 ) ,P(5) 

IN= ( I + 1 ) *N 
JN=(J+1)*n 
R( I t J*3)=P( 5) 

R(I»J»4) = (P(1 ) -P ( 2 ) )/L( ID 
R ( I,JtD=P(2)+JN*R( I,J,4) 

R ( I »J»2 )=P(3)+(P(A)-IN— 1 )*P( 5) 

R ( I t J » 5 ) =L ( 29 ) 

RETURN 

END 


RAYSOOO 1 
RAYS0002 
RAYSOOO 3 
RAYS0004 
RAYS0003 
RAYS0006 
RAYS0007 
RAYS0008 
RAYS0009 
RAYSOOIO 
RAYSOOl I 



SUBROUTINE RDPIK ( I NA, INB, KOA , KOB ,’L' ,'N','R , I A', JA, I B, J8, IK,JK, IR»JRtP 
ltJDISP) 

DIMENSION INA( IA,JA),INB( IB,JB),KOA( IK,JK),KOB(lK,JK) 

DIMENSION L(32 ) ,R< IR, JR,5) ,P(5 ) 

DO 400 J=l, JR 
DO 300 1=1 ,IR . 

CALL RAYSET ( L , R , I , J , P , N , I R., JR ) 

WRITE (15) (R( I ,J,K ) ,K=1,5> 

JDL=( J-l )*N+JDISP 
I DL= ( 1-1 ) *N 
DO 200 JX=1,JK 
J Y= JDL+ JX 
DO 100 I X = 1 » IK 
I Y= I DL+ 1 X 

IF ( IX.GT.U 10) ) GO TO 20 

IF ( JX.GT.L( 11 ) ) GO TO 20 • * 

10 CONTINUE 

KOA( IX,JX)=INA( IY,JY) 

KOB( IX, JX )=INB( I Y, JY) 

GO TO 30 
CONTINUE 
OA ( IX,JX)=0 

KbB( IX, JX )=0 \ 

CONTINUE '\^ 



30 

100 

200 


CONTINUE 
CONTINUE 
WRITE (15) KOA 
WRITE ( 15) KOB 
300 CONTINUE 
400 CONTINUE 

READ (13J.INA 
READ (14) INB 
RETURN 
END 


RDPIOOO. 

rdpiooo; 

RDPIOOO. 
RDPIOOO' 
RDPIOOO: 
RDPIOOOt 
RDPIOOO 
RDPIOOO.' 
RDPIOOO' 
RDPIOOK 
RDPI001 
RDPI001 
RDPI001 
RDPI001 
RDPI001 
RDP 1 0 0 1 < 
RDP 1 001 
RDP 1 00 1 : 
RDPIOOl' 
RDPI002* 
RDP 1 002 
RDP 1002 
RDP 1 002 
RDP 1002* 
RDP 1 002: 
RDP 1002 
RDP 1 002 
RDP 1002 
RDP 1 002 
RDP 1003: 
RDP 1 003 
RDP 1003 
RDP 1 003 
RDP 1 003 
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SUBROUTINE RDREC (N,N,I,J) 
DIMENSION M( I ) ,N< J ) 

READ < 13 ) M 
READ (14) N 
RETURN 
END 


RDREOOO 

RDREOOO 

RDREOOO 

RDREOOO 

RDREOOO 

RDREOOO 
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SUBROUTINE RDTAPE RDTAOOOl 
C SUBROUT INE FOR READING TAPES FOR . SUCCESS I VE ORBITS AND SAME ARE ARDTA0002 
C AND FORMING PAIRS OF 64X64 ARRAYS. RESULTS ARE OUTPUT ON A SINGLE RDTAOOOB 
f TAPE IN THE SEQUENCE OF ONE RECORD CONTAINING TITLE, LOCATION AND RDTA0004 
(. NUMB.ER OF ARRAY PAIRS FOLLOWED BY TRIPLES OF ONE RECORD WITH TIT LERDTAOOO 5 
C I,J, AND APPROXIMATE COORDINATES, ONE RECORD OF DATA FOR THE 64X64 RDTA0006 
C - POINTS FOR THE- FIRST PASS, AND QNE 'RECORD OF DATA FOR THE SECOND RDTA0007 
C PASS. THE LAST TRIPLE IS FOLLOWED BY A RECORD OF THE TITLE, PLACE RDTA0008 
(, NUMBER OF PAIRS AND THE NMR TAPE DDNAME, STARTING AND ENDING T I ME SRDTA0009 
C 


AND OTHER DATA 


FOR EACH OF THE TAPES 
COMMON INTEGRI 88192 ) 

DIMENSION INDATA (40000) , I NDATB ( 40000 ) 

KOUTA(4096) ,K0UTB(4096) 

INBUFAI 26 ) , BUF A ( 7 ) ,P ICTAI 5 ) , NGR I DA ( 4 ) , TNMRA ( 8 ) , XNMRA ( 7 ) 
INBUFBI 26) ,BUFB( 7) ,PICTB( 5) , NGR I DB ( 4 ) , TNMRB ( 8),XNMRB(7) 
AB F ( 18) ,BBF( 18) , L I ST ( 3 2 ) , T I TLE ( 3 ) , PL ACE ( 5 ) , NPRS ( 3’) 

PARM ( 20) ,RAYS( 500) 

( I ND AT A ( 1 ) , INTEGRI 1 ) ) , ( I NDATB ( 1 > , I NTEGRI 4000-1 ) ) 

I KOUTAI 1 ) , INTEGRI 80001 ) ) , I KOUT B I 1 ) , I NTEGR I 84097 ) ) 

I INBUFAl 1),BUFA( 1) ),( INBUFAI8) , XNDRA ) , I INBUFAI 9) 

I INBUFAI 10 ) ,P ICTAI 1 ) ) , I INBUFAI 15) , NGR I DA I 1 ) ) 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
UI VALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 


E9 


RDTA0010 
RDTA001 1 
RDTAOO 1 2 
RDTAOO 1 2 
RDTAOO 14 
RDTA001 5 
RDTAOO 16 
RDTAOO 1 7 
RDTAOO IE 
RDTAOO 1 9 
,MRCA )RDTA002C 
RDTA0027 


I INBUFBI 1),BUFB( 1) ),( INBUFBI8 ) , XNDRB ) , I INBUFBI 9) ,MRCB ) RDTAOO 22 


RDTA0022 
RDTA0024 
RDTA002 5 


I INBUFBI 10 ) ,PICTB( 1 ) ) , I INBUFBI 15) ,NGRIDB< 1 ) ) 

I INBUFAI 19) , TNMRA I 1 ) ) , I INBUFBI 19) , TNMRB I 1 ) ) 
(LIST(l),KEY),(eiST(2),TITLE(l) ) , I L I ST I 5 ) , PL ACE 1 1 ) ) 

I LIST! 10) tNPRSI 1 ) ) , I LIST! 13) , XNMRA I 1 ) ) , I L I ST I 32 ) , KEND ) RDT A0026 
I LISTI20),XNMRB( 1) ) , I L I ST I 27 ) , XNAD I R ) , I L I ST I 28 ) , MERC ) RDTA002 7 
(L 1ST! 29 ) ,NTDIF ) , I INBUFAI 1 ) , ABF I 1 ) ) , I INBUFBI 1) ,BBF( 1 ) ) RDTA002E 
(LI ST I 30) , IRAYS ) » (L 1ST I 31 ) »JRAYS) 


EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 

1 FORMAT I 20A4 ) 

2 FORMAT (75H1RUN QUESTIONED 
1MERCAT0R PROJECTION. ) 

3 FORMAT (71H1RUN QUESTIONED 
IE NOT IDENTICAL. ) 

4 FORMAT I 65H1RUN TERMINATED 
1TA POINTS. ) 

5 FORMAT I 46H0 THE TIME DIFFERENCE 

1MINUTES. ) 

6 FORMAT I 36H0 ARRAY ORGANIZATION STEP 

7 FORMAT I 16X, 1 4 , 6 0 X ■) 

8 FORMAT I 1 HI , 20A4 ) 

9 FORMAT I 16H0NUMBER OF ROWS= ,16, 20H 
100 CONTINUE 

REWIND 13 
REWIND 14 


BECAUSE AT LEAST ONE INPUT TAPE WAS 


BECAUSE PICTURE LOCATIONS OR SCALES 


BECAUSE INPUT AREA LARGER THAN 


RDTA002 C . 
RDTA003C 
N0TRDTA003) 
RDTA0032 
WERRDTA0032 
RDT A003^ 
40000 DARDTA003 5 
RDT A003( 


FOR THIS PICTURE PAIR IS,I8,10H 


COMPLETED. 


NUMBER OF COLUMNS= ,16) 


RDTA003*. 
RDTA003E 
RDTA0039 
RDTA004C 
RDTA004 ) 
RDTA004; 
RDTA004I 
RDTA004^ 
RDTA004I 
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ND I SP=1 6 

READ (5,1) PARM 

READ (5,7) NT I ME 

PRINT-8, P ARM - 

DO 110 L=1 ,3 
T I TLE ( L ) =PARM ( L ) 

110 CONTINUE 

READ (13) ABF 
READ (14) BBF 

IF ( MRCA. EQ.MRCB ) GO TO 120 
PRINT 2 
120 CONTINUE 

DO 130 LA = 1 , 5 

IF (PICTA(LA).NE.PICTB(LA) ) GO TO 150 
130 CONTINUE 

DO 140 LB = 2, 3 

IF (NGRIDA(LB).NE.NGRIDB(LB) ) GO TO 150 
140 CONTINUE 

C PROJECTIONS, LOCATIONS AND SCALES OF TWO SAMPLES MATCH. 
GO TO 200 
150 PRINT 3 
200 CONTINUE 
KE Y = 1 

DO 210 LC = 1 , 5 
PLACE(LC)=PICTA(LC) 

210 CONTINUE 

NPRS( 1 )=NGRIDA(2 ) 

NPRS( 2 ) = NGRIDA( 3 ) 

NPRS ( 3 ) =NPRS ( 1 ) ^NPRS < 2 ) 

PRINT 9»NPRS(2),NPRS( 1 ) 

XNAD I R=XNORA 
MERC=MRC A 

IRAYS=(NPRS( 1 )/NDISP)-3 

JRAYS=( NPRS(2 J/NDISP )— 3 

IF ( IRAYS.LT. 1 ) I R A Y S = 1 

I F ( JRAYS.LT. 1 ) JRAYS= 1 

IF (NPRS(3 ) .LE. 40000) GO TO 220 

PRINT 4 

STOP 

220 CONTINUE 

JD I SP =N PRS ( 2 )-{JRAYS + 3)*NDISP 
KEND= 1 

WRITE (15) LIST 
I A = L I ST ( 10) 

JA = L I ST ( 1 1 ) 


RDTA0046 
RDTA004 7 
RDTA004E 
R D T A 0 0 4 c . 
RDTA005C 
RDTA005 1 
RDTA0052 
RDTA005 2 
RDT A0054 
RDTA0052 
RDTA005L 
RDTA005 "i 
RDTA005E 
RDTA005S 
. RDTA006C 

RDTA006 ] 
RDTA0062 
RDTA0062 
FORM MATCHINGRDTA006' 
RDTA006I 
RDTA0066 
RDTA006T 
RDTA006E 
RDTA006L 
RDTA007C 
RDTA007 ] 
RDTA0072 
RDTA0072 
RDTA007' 
RDTA007 i 
R0TA007( 
RDTA007' 
RDTA007I 
RDTA007 1 
RDTA008C 
RDTA008 I 
RDTA00 82 
RDTA008 
RDTA008' 
R D T A 0 0 8 ‘ 
RDT A008< 
RDTA008 ' 
RDTA008! 
RDTA008 ‘ 
RDT A009( 
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I B = L I S.T ( 10) 
JB=LIST.( 11 ) 
IK=NDI SP*4 
JK = ND I S P-4 
I R = L I ST ( 30) 
JR = L I ST (31 ) 
I RD= I A*J A 
JRD= I B* JB 


RDTA009 
RDTA00U, 
RDTA009 . 
RDT A009' 
RDTAOOV 
RDT A00y< 
RDTA009 ' 
RDTA009? 


CALL RDREC ■ < I NDATA, I NDATB , I RD , JRD ) RDTA009' 

300 CONTINUE RDTA010( 

310 CONTINUE RDTAOIO: 

NTD I F =NT I ME ’ RDT AO 1 01 

PRINT 5 , NTD I F RDTAOIO: 

CALL RDPIK ( INDATA, INDATB ,KOUT A , KOUTB , L I ST ,ND I SP , RAYS , 1 A, JA, I B , JB , RDT AO 1 0< 
1 1 K , JK , IR, JR, PLACE, JDISP) RDT AO 10 i 

REWIND 13 * RDT AO 1 Of. 

REWIND 14 RDTAOIO" 

L I ST ( 1 ) =9999 ’ RDTAOIO! 

WRITE (15) LIST - RDTA010 C 

END FILE 15 ' RDTA011C 

REWIND 15 -RDTA01 1 1 

,RINT 6 RDTAOli; 

5000 CONTINUE v RDTA011- 

RETURN \ RDTA0114 

END \ \ RDTA0115 
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SUBROUTINE SETHKM ( NHARM , I FE RR , I F S , M , MM, NN, U, I NV, S ) 

SETHOOO 1 


DIMENSION U(MM,NN ) ,M(3 ) , I N V ( MM ) , SIMM) 

SETHOOO; 


COMPLEX U 

SETHOOOl- 


DO 600 I =1 ,MM 

SETHOOO*- 


SC I ) = o 

SETHOOO f 


I NV ( I ) =0 J 

SETHOOOL 

600 

CONTINUE 

SETHOOO' 


CALL HARM (U,M,INV,S»0»IFERR) 

SETHOOO/ 


IF ( I F ERR ) 601,602,601 

SETHOOO" 

601 

CALL IFEXIT (0 , I FERR , I F S , M , MM , NN , U , I NV , S ) 

SET HO 0 1 f 

602 

CONTINUE 

SETH001 i 

# 

RETURN 

SETHOO 1 2 


END 

SETHOOll 


o o o o 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE SETRAY ( RL , DEL I , DE L J * T I ME K ) 


I KUUlJ J 


INPUT TO THIS SUBROUTINE IS THE FIVE-WORD RECORD READ FROM ARRAY 
TAPE PRECEDING THIS ARRAY PAIR. IT CONTAINS ARRAY-PAIR VALUES OF 
LATITUDE IN DEGREES NORTH LATITUDE 
LONGITUDE IN DEGREES WEST LONGITUDE 
LATITUDINAL SCALE IN DEGREES 
LONGITUDINAL SCALE IN DEGREES 

TIME BETWEEN DATA MEASUREMENTS OF WINDOW AND SEARCH DATA 
(TIME IS IN MINUTES) 

OUTPUT FROM THE SUBROUTINE ARE THE THREE CONVERSION PARAMETERS 
NEEDED TO CONVERT THE LAG COORDINATES OF THE CROSS-CORRELATION 
PEAK INTO A MOTION VECTOR FOR THE WINDOW USED 

DELI = THE EAST-WEST SPACING UNITS IN NAUTICAL MILES 
DEL J = THE NORTH-SOUTH SPACING UNITS IN NAUTICAL MILES 
T I MEK= THE TIME BETWEEN DATA MEASUREMENTS IN HOURS 


DIMENSION RL( 5) 
RADS=RL ( 1 ) / 57 .295 
DEL I = ( R L ( 3 ) * 60.0) 
DEL J = RL ( 4 ) * 60.0 
T I MEK= RL ( 5 ) /60.0 
RETURN 
END 


COS( RADS ) 
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•SUBROUTINE SPLITV ( MM , NN , MMNN , K I I , UU t V V ) 
DIMENSION UU(MMNN) ,VV(MMNN) 

COMPLEX UU, VV,UP,UM,UX,UY 
UU( 1 ) =REAl ( V V ( 1 ) )*AIMAG( VV( 1 ) ) 

I I I = MM + 2 
DO 3001 1=2, MM 
I I = I I I - I 
UP=VV( I ) 

UM= V V ( II) 

UX=CONJG(UP )+UM 
UY=UP— CONJG(UM) 

UU( I ) = (0.,-0.25)*UX*UY 

3001 CONTINUE 
MMNN=MM*NN 
JB=MM+ 1 
JE=MMNN-MM+1 
JJJ=MMNN+2 

DO -3002 J = JB , J E , MM 

JJ=JJJ-J 

UP = VV ( J ) 

UM= VV ( J J ) 

UX=CONJG( UP ) +UM 
UY=UP-CONJG(UM) 

UU( J)=(0.,-0.25)*UX*UY 

3002 CONTINUE 
JB= JB+1 
JE= JE+ 1 
MMM=K I 1-1 

I I I =MMNN+ I I I 
DO 3004 J = JB , JE ,MM 

I E = MMM+ J 

DO 3003 I -J , I E 

I I = I I I - I 
UP = VV( I ) 

UM = VV ( II) 

UX=CONJG(UP )+UM 
UY=UP-CONJG(UM) 

UU( I)=(0.,-0.25)*UX*UY 

3003 CONTINUE 

3004 CONTINUE 
JB = J B + K I I 
JE = JE+K 1 1 
MMM=MMM- 1 

DO 3006 J = JB , JE ,MM 

I E=MMM+ J ' _ 


SPL 1000) 
SPL 1 0002 
SPL I 0001 
SPL 1000* 
SPLI0001 
SPL 1 0006 
SPL I 0001 
SPLIOOOt 
SPLIOOO 1 . 
SPL IOOK 
SPLI001) 
SPLI0011 

spliooi: 

SPL 1001* 
SPL I 001 1 
SPL 10016 
SPL 1001 ~ 
SPLIOOH 
SPL 1 00 1 ‘ 
SPL 100 2C 
SPL 1 002) 
SPLI002; 
SPLI002: 
SPL 1002* 
SPLI002! 
SPL 1 002< 
SPL 1002 
SPL 1002) 
SPL I 002‘ 
SPLI003C 
SPLI003: 
SPL 1003. 
SPLI003 
SPL 1003- 
SPL I 003: 
SPLI003* 
SPL 1003 
SPLI003) 
SPLI003* 
SPL 1 004( 
SPL I 004 
SPL 1004. 
SPL 1 004 
SPL 1004- 
SPL I 004 


DID 3005 I = J » I E SPL 10046 

I I = I I I - I SPL I 004 7 

UU( I ) =C0NJG( UU ( 1 I ) ) SPL 1 0048 

3005 CONTINUE $PLI0049 

3006 CONTINUE SPLI0050 

RETURN SPL I 005 1 

END-;, ' SPL 10052 


> 




SUBROUTINE THRESH ( M ARR A Y ,MM ,NN , M I NCL D , MAXCLD , NFREO , PRC ENT , XMEAN , THREOOOl 
1SDEV, TARRAY, INIT I X,NDEV ) THRE0002 

DIMENSION TARRAY(MM,NN) ,M ARRAY I MM , NN )■ THRE0003 

I FORMAT I 20H0THRESH0LDS ARE MIN=,I4,5H MAX=,I4,21H VALID POINTS NUMTHREOOQ4 



1BER * I 7 , 2 1 H FOR A PERCENTAGE OF ,F6.2,10H THE ME AN= , F 7 . 1 , 20H 

STATHRE0005 


2NDARD DEVIATI0N=,F7.1) 

THRE0006 


M S U M = 0 

THRE0007 


MSQSUM=0 

THRE0003 


NFRE0=0 

THRE0009 


PKCENT = 0 

THREOO 1 0 


DO 20 J=1,NN ■ 

THRE0011 


DO 10 1=1, MM 

THREOO 1 2 


IF (MARRAYI I ,J ) .LT.MINCLD) GO TO 10 

THREOO 1 3 


IF (MARRAYI I, J).GT. MAXCLD) GO TO 10 

THRE001 4 


NFREQ=NFREQ+1 

THREOO 1 5 


MSUM = MSUM + MARRAYI I , J ) 

THREOO 1 6 


MSOSUM=MSOSUM + MARRAYI I , J ) * * 2 

THREOO 17 

10 

CONTINUE * 

THRE001 8 

20 

CONTINUE 

THREOO 1 9 


XFREO=NFREO 

THRE002 0 


_XSUM=MSUM 

THRE002 1 


XSQSUM=MSQSUM 

THRE002 2 


XMEAN = XSUM/XFREQ 

THRE0023 


PRCENT = (XFRE0*100.0)/IMM-NN) 

THRE0024 


XVAR= I ( XSOSUM/XFREO)-! XMEAN=**2 ) ) *( XFREQ/ I XFREQ-1 .0 ) ) 

THRE0025 


SDEV = SORT I XVAR ) 

THRE002 6 


DO 60 J=1,NN 

THRE0027 


DO 50 1=1, MM 

THRE002 8 


IF (MARRAYI I , J ) . LT . M I NCLD ) GO TO 30 

THRE0029 


IF (MARRAYI I ,J).GT. MAXCLD) GO TO 31 

THRE0030 


TARRAY I I , J ) = MARRAYI I, J) 

THRE0031 


GO TO 40 

THRE0032 

30 

CONTINUE 

THRE0033 


DEVN=NDEV/( 6.0 ) 

THRE0034 


ADEV=SDEV*DEVN 

THRE0035 


GO TO 32 

THRE0036 

31 

ADE V=0 • •- 

THRE0037 

32 

CALL BKGRND I I N I T I X, ADE V, XME AN, V ) 

THRE0038 


TARRAY ( I , J ) =V 

THRE0039 

40 

CONTINUE 

THRE0040 

50 

CONTINUE 

THRE0041 

60 

CONTINUE 

THRE0042 


NSUBS=(MM*NN )-NFREQ 

THRE0043 


PRINT 1,MINCLD, MAXCLD, NFREQ, PRCENT , XME AN, SDEV 

THRE0044 


RETURN THRE0045 
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END 


THRE0046 
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MEANS=»F8.3,48H 
MEANS=,F8.3 > 

STANDARD DE V I AT I ONS= , F8 . 3 , A8H 
DEVIATI0NS=»F8.3 ) 



SUBROUTINE T 1MNSD ( Y, YWMN , YWSD , MM, NN , MW , NW, YMNMN, YMNSD, 
1ZTEST) 

DIMENSION YWMN(MW,NW) , YWSD(MW,NW) ,Y(MM,NN) 

MUST HAVE MW= ( MM/2 ) +1 » NW=(NN/2)+l UPON CALLING 

1 FORMAT ( 123HOFOR THE SEARCH AREA SUBARRAYS THE FOLLOWI 

IS WERE FOUND AFTER THE SEARCH AREA .MEAN WAS SUBTRACTED 
2TA ) 

2 FORMAT ( 3 1 HO MEAN OF 
1DEVIATI ON OF 

3 FORMAT ( 3 1 HO MEAN OF 
1 DEV I AT I ON OF STANDARD 

ZTEST=1 
YMNMN=0 
YMNSD=0 
YSDMN=0 
YSDSD=0 
YSUM=0 
YSQSUM=0 
MWN=MW*NW 
M=MM/2 
NN/2 
KN=M*N 
z='mn 
ZW = M'WN 

zs=z*z 

ZWS=ZW*ZW 
DO 10 1=1, M 
DO 10 J = 1 , N 
YSUM=YSUM+Y ( I , J ) 

YSQSUM=YSQSUM+Y( I , J )#Y( I , J) 

10 CONTINUE 
■ YWMN ( 1 1 1 ) =YSUM 
YWSD (1,1 ) =YSQSUM 
SQUARE IN LOWER LEFT QUADRANT 
NOW PROCEED UP THE 
DO 30 J=1,N 
JM= J 
JP= J+N 
JW=J+1 

YSUM=YWMN { 1 , J ) 

YSQSUM=Y WSD (1,J) 

DO 20 1=1, M 

YSUM= YSUM+Y (I,JP)-Y(I,JM) 

YSOSUM=YSQSUM+Y( I , JP ) #Y ( I , JP ) -Y ( I , JM ) *Y ( I , JM ) 

20 CONTINUE 




\ 


PROCESSED AND STORED 
COLUMN ALONG THE LEFT SIDE 


YSDMN, YSDSD,T1MN0001 
T 1MN0002 
T1MN0003 
T1MNOOOA 
NG STATI STICT1MN0005 
FROM THE DAT 1 MN0006 
T 1MN0007 
STANDARD’ T1MN0008 
T1MN0009 
STANDARD T1MN0010 
T1MN001 1 
T1MN0012 
T1MN0012 
T1MN001A 
T1MN0015 
TIMNOOlfc 
' T1MN0017 
T1MN0018 
T1MN001S 
T1MN0020 
T1MN002 1 
T1MN0022 
T1MN0022 
T 1MN002A 
T1MN002 5 
T1MN0026 
T1MN0027 
T1MN0028 
T1MN002? 
T1MN0030 
T1MN003 1 
T1MN0032 
T1MN0032 
T1MN003A 
T1MN0035 
T1MN0036 
T1MN0037 
T1MN0036 
T1MN0039 
T1MN00AC 
T1MN00A) 
T1MN00A2 
TIMNOOAr 
TIMNOOA^ 
T1MN00AS 
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YWMN ( 1 t JW ) = YSUM 
YWSD ( 1 , JW)=YSOSUM 
30 CONTINUE 

C PROCESSING LEFTMOST COLUMN OF M ROWS COMPLETED 
C CONTINUE PROCESSING SHIFTING ONE COLUMN TO RIGHT AT A TIME 
CC AND GETTING RESULTS FOR* A NE W STR I P M ■ W I DE BY N ROWS EACH TIME 
DO 70 1STEP=1,M 
I M= I STEP 
I P= I STE P+M 
IW=ISTEP+1 
YSUM=YWMN ( I ST E P » 1 ) 

YSOSUM=YWSD (ISTEP,1) 

DO 40 J =1 »N 

YSUM=YSUM+Y( I P, J ) — Y ( IM, J ) 

YSQSUM=YSQSUM+Y( IPr J)*Y( IP,J)-Y( IM, J)*Y(IM,J) 

40 CONTINUE 

YWMN ( I W 1 1 ) = YSUM . 

YWSD( IWt 1 )=YSQSUM - * . 

C BOTTOM M BY N BOX OF COLUMN PROCESSED 
C NOW COMPLETE PROCESSING REMAINDER OF COLUMN 
DO 60 J=1,N 
JM= J 
JP= J+N 
JW=J+1 

YSUM=YWMN( IW T J ) 

YSOSUM=YWSD( I W t J ) 

I LO=I STEP+1 
I HI= I STEP+M 
DO 50 I =ILOt IHI 
YSUM=YSUM+Y ( I , J P ) -Y ( I , JM ) 

YSQSUM = YSQSUM + Y ( I » JP ) *Y ( I » JP ) — Y ( I » JM ) *Y ( I » JM ) 

50 CONTINUE 

YWMN (IW t JW)=YSUM 
YWSD ( IW t JW)=YSOSUM 
60 CONTINUE 

C THIS COLUMN HAS BEEN PROCESSED. GO BACK TO START A NEW ONE. 

C COLUMN INCREMENTED TO RIGHT BY ONE ELEMENT-COLUMN 
70 CONTINUE 

C ALL (M+l) BY (N+l) M BY N ARRAYS ARRAYS PROCESSED . . 

C ALL ( M+l )X ( N+l ) M BY N ARRAYS PROCESSED TO THE EXTENT OF GETTING 
C SUMS AND SUMS OF SQUARES. CONVERT THESE TO MEANS AND SIGMAS 
DO 80 I = 1 t MW 
DO 80 J = 1 * NW 
XS = Z*YWSD( I » J ) 

SXS = YWMN ( I ,J )*YWMN( I , J) 


T1MN0046 
T1MN0047 
T1MN0048 
T1MN0049 
T1MN0050 
T1MN0051 
T 1MN0052 
T1MN0053 
T 1MN0054 
T1MN0055 
T 1MN0056 
T1MN0057 
T1MN0058 
T 1MN0059 
T1MN0060 
TIMN006I 
T1MN0062 
T 1MN0063 
T 1MN0064 
T1MN0065 
T 1MN0066 
T1MN0067 
T1MN0068 
T1MN0069 
T1MN0070 
T1MN0071 
T1MN0072 
T1MN0073 
T1MN0074 
T1MN0075 
T1MN0076 
T1MN0077 
T 1MN0078 
T1MN0079 
T 1MN0080 
T1MN0081 
T 1MN0082 
T1MN0083 
T 1MN0084 
T1MN0085 
T1MN0086 
T 1 MN008 ’/ 
T1MN0088 
T1MN008S 
T 1MN009C 
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YWV=(XS-SXS)/(ZS-Z) 
YWM=YWMN ( I » J ) /Z 
YMNMN= YMNMN+ Y WM 
YWMN ( I » J ) = YWM 
YMNSD=YMNSD+( YWM*YWM ) 
IF( YWV.LE.O) GO TO 75 
YWSD ( I ,J)=SQRT( YWV) 
YSDMN=YSDMN+YWSD( I , J ) 
YS[)SD=YSDSD + YWV 
GO TO 80 

75 YWSD ( I , J ) =— 1.0 
ZTEST =0 

80 CONTINUE 

XMM=YMNMN/ZW 
XMN=YMNMN*YMNMN 
XMV=YMNSD*ZW-XMN 
XMV=XMV/( ZWS-ZW) 

IF (XMV.LE.O) GO TO 85 
XMS= SORT ( XMV ) 

- GO TO 86 

85\XMS=- 1 • 0 

86 YNNSD=XMS 
YMNMN=XMM 
XSM=^SDMN/ZW 
XSN= YSDMN^YSDMN 
XSV=YSDSD#ZW— XSN 
XSV=XSV/( ZWS-ZW) 

IF (XSV.LE.O) GO TO 87 
XSS= SORT ( XS V ) 

GO TO 88 

87 XSS=-1.0 

88 YSDSD=XSS 
YSDMN= XSM 
PRINT 1 

PRINT 2» YMNMN, YMNSD 
PRINT 3 t YSDMN » YSDSD 
RETURN 
END 


T1MN0091 
T1MN0092 
T1MN0092 
T1MN009* 
T1MN0095 
T1MN0096 
T 1MN009 / 
TlMN009f 
T 1 MN009 C . 
T1MN010C 
T1MN010) 
T1MN010Z 
T1MN010I- 
T1MN010 Z 
T1MN0101 
TIMNOlOf 
TIMNOIOI 
T1MN010E 
T1MN010 4 
TlMNOllt 
T1MN011 3 

timnou; 

timnou: 

T1MN011* 
TIMNOllf 
T1MN01H 
T1MN011'/ 
TlMNOllt 
TIMNOll 4 
T1MN012C 
T1MN012) 
T1MN012; 
T1MN012; 
T1MN012' 
T1MN0I2! 
T 1MN0 1 2< 
T1MN01 2 ' 
T1MN0121 



SUBROUTINE VECTOR ( I BM , JBM , DEL I * DE L J 
7 FORMAT (23HO THE PEAK VALUE AT 1= , I 
IX MOTION WAS FROM ,I3,12H DEGREES AT 
2 I = I BM 
Z J= JBM 
W I = Z I *DE L I 
WJ=ZJ*DEL J 
WSQ=WI*WI+WJ*WJ 
WVEL=SORT(WSO)/TIMEK 
MVEL = WVE L 
.IF (JBM) 44,40,44 

40 *IF ( IBM ) 41,42,43 

41 MPH I =90 
GO TO 70 

42 MPH 1=0 
GO TO 70 

43 MPH 1=270 
GO TO 70 

M WPH I =W I / W J 

VAPHI=ABS< WPHI ) 

BPHI=ATAN( APHI ) 

PHI=BPHI*( 57.295) 

NPH I = PH I 

IF (JBM) 45,51,51 

45 IF (IBM) 46,47,48 

46 MPHI = NPH I 
GO TO 70 

47 MPHI =360 
GO TO 70 

48 MPHI=360-NPHI 
GO TO 70 

51 IF (IBM) 52,53,54 

52 MPH 1=180 — NPHI 
GO TO 70 

53 MPH 1=180 

GO TO 70 „ 

54 MPH I = 1 80+NPH I 
70 CONTINUE 

MPK=AMAX*400.0 

PRINT 7, IBM, JBM, MPK, MPHI, MVEL 

RETURN 

END 


T I MEK , AMAX ) VECT0001 

,9H AND J= , I 3 , 4H IS, 19, 18HVECT000 2 
,14, 9H KNOTS ) VECT0003 

VECT 0004 
VECT0005 
VECT0006 
VECT0007 
VECTOOO 8 
VECT0009 
VECT0010 
VECT001 1 
VECT001 2 
VECT0013 
VECT0014 
VECT 00 1 5 
VECT0016 
’ VECTOO 1 7 

‘ VECTOO 1 8 

VECTOO 1 9 
VECT0020 
VECT0021 
VECT002 2 
VECT0023 
VECT 0024 
VECT0025 
VECT 002 6 
VECT0027 
VECT002 8 
VECT0029 
VECT0030 
VECT0031 
VECT0032 
VECT0033 
VECT0034 
VECT0035 
VECT0036 
VECT0037 
VECT003 8 
VECT0039 
VECT0040 
VECT0041 
VECT 0042 



ooooooooooo 


SUBROUTINE WINDOW 
1NW,SDEVA,MXSWCH) 


INPUTS ARE 


( WT,MWT,NWT,XCV,MVT,NVT, XC VP , MPT , NPT , YSDLAG,MW, 


WT, INVERSE TRANSFORM OF W IN LAG FORM AND ITS DIMENSIONS 
YSDLAG, THE STANDARD DEVIATIONS OF SEARCH' AREA SUBARRAYS 
SDEVA, THE STANDARD DEVIATION OF THE WINDOW SUBARRAY 
WORK AREA XCV AND OUTPUT ARRAY AREA AND THEIR DIMENSIONS 

OUTPUT WILL BE A CROSS-CORRELATION ARRAY XCVP FOR LAGS OF THE WINDOW 


TO CALL THIS SUBROUTINE 


SET MVT = ( MWT-1 ) /2» NVT = ( NWT-1 ) /2 
AND MPT = MVT+1 NPT = NVT+1 


C 

C 


C 

C 

C 

C 

c 

c 


DIMENSION WT(MWT,NWT), XCV< MVT ,NVT ) , XCVP ( MPT » NPT ) » YSDLAGI MW ,NW ) 
COMPLEX WT 

1 FORMAT ( 83HOCROSS CORRELATIONS RECOMPUTED USING UNIFORM STANDARD 
1 EV I AT I ON OVER SEARCH AREA. ) 

I NC-MVT J 2 
JNC=NVT / 2 

FORM XCV BY SELECTING REAL PARTS OF WT FOR CENTRAL FOURTH ONLY 
DO 100 1=1, MVT 
DO 100 J=1,NVT 
I W= I + 1 NC 

JW = J + JNC - ■■ 

XCV( I ,J )=REAL(WT( IW,JW) ) 

XCVP( I , J )=XCV( I , J ) 

100 CONTINUE 

COMPUTE CROSS-CORRELATION COEFFICIENT S FOR NORMALIZED SUBARRAYS BY 
DIVIDING CROSS-COVARIANCE FOR EACH LAG PAIR BY STANDARD DEVIATION 
OF WINDOW SUBARRAY AND THE STANDARD DEVIATION OF THE SEARCH AREA 
SUBARRAY FOR THAT PAIR OF LAGS 

IF { MXSWCH.NE.O ) GO TO 500 
DO 400 1=1 ,MPT 
DO 400 J = 1 , NPT 
DEN=SDEVA*YSDLAG( I , J ) 

IF ( DEN. LT.0.01 ) GO TO 500 
400 XCVP( I,J)=XCVP(I,J)/DEN 
RETURN 

500 DEN=SDEVA*SDEVB 

IF ( DEN. LT.0.01 ) GO TO 700 
DO 600 1=1, MPT 
DO 600 J=1,NPT 


WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
DWI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
. WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 
WI 


ND0001 
ND0002 
ND0003 
ND0004 
ND0005 
ND0006 
ND0007 
ND0008 
ND0009 
ND0010 
ND001 1 
ND0012 
ND0013 
ND0014 
ND0015 
ND0016 
ND0017 
ND0018 
ND0019 
ND0020 
ND0021 
ND0022 
ND0023 
ND0024 
ND0025 
ND0026 
ND0027 
ND0028 
ND0029 
ND0030 
ND0031 
ND0032 
ND0033 
ND0034 
ND0035 
ND0036 
ND0037 
ND0038 
ND0039 
ND0040 
ND0041 
ND0042 
ND0043 
ND0044 
ND0045 



600 

XCVP( I ,J )=XCVP(I »J)/DEN 

WIND00A6 


PRINT 1 

WIND00A7 


RETURN 

WIND00A8 

700 

PRINT 1 

WIND00A9 


PRINT 2 

WIND0050 


RETURN 

WIND0051 

2 

FORMAT ( 53H0C0MPUTAT ION REJECTED BECAUSE CLOUD DATA TOO SPARSE 

)W IND0052 


END 

W1ND0053 
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( IH ,30X, I 5 1 3X » I5*8X, 15) 

( 132H0THIS. MOT I ON . COMPUTAT ION WAS REJECTED 
LAG PAIRS HAD GREATER THAN 0.99 CROSS CORR 

) 

( 1 05H0THI S MOTION COMPUTATION WAS REJECTED 
CORRELATION COEFFICIENT EXCEEDED 1.005. 


SUBROUTINE XCMAX ( A , AM A X , MM , NN , I MAX , JMAX , MXS WCH ) 
DIMENSION A ( MM » NN ) » VM AX( 6 ) tMAX I ( 6 ) » MAX J ( 6 ) t MAX V ( 6 

1 FORMAT (58HOEXTREME MAXIMA NOTED WERE FOR I LAG = 

IE ) 

2 FORMAT 

3 FORMAT 
IN FIVE 

. 2NT . 

4 FORMAT 
1 CROSS 

NMAX=0 
MXSWCH=0 
ATST = 0. 247 5 
ALIMIT=0. 25125 
DO 10 K= 1 * 6 
VMAX ( K ) =0 
MAXI (K 1=0 
MAX J ( K ) =0 
MAX V ( K ) = 0 
.CONTINUE 
JX=0 


XCMAOOO 1 

) XCMA0002 

J L AG= THE VALUXCMAOOO 3 

XCMA0004 

XCMA0005 

BECAUSE MORE THAXCMA0006 
ELATION C0EFFICIEXCMA0007 

XCMAOOO 8 

BECAUSE C0MPUTEDXCMA0009 


) 



11 

12 


128 

129 

130 


AM = t 
BM=0 N 


x. 


\ 


DO 23 1=1 *MM 

DO 23 J=1,NN 

IF ( A( I, J) ) 13*13*11 

IF ( A ( I , J ) - AM ) 13*13*12 

AM= A ( I,J ) 

1 MX= I 
JMX= J 

IF (AM.LE.ATST) GO TO 130 

IF ( AM.GE.AL IMIT ) GO TO 128 

NMAX=NMAX+1 

VMAX ( NMAX ) =AM 

MAXI ( NMAX ) = I — 1 — < MM/2 ) 

MAXJtNMAX )=J— l-INN/2 ) 

MAXVINMAX )=VMAX( NMAX) *400.00 

IF (NMAX. IE. 5) GO TO 130 • 

PRINT 3 

GO TO 129 

PRINT 4 

MXSWCH=1 

GO TO 100 

CONTINUE 
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XCMAOO 10 
XCMA001 1 
XCMAOO 12 
XCMA001 3 
XCMA0014 
XCMAOO 1 5 
XCMAOO 1 6 
XCMA001 7 
XCMAOO 1 8 
XCMAOO 1 9 
XCMA0020 
XCMA002 1 
XCMA0022 
XCMA002 3 
XCMA0024 
XCMA002 5 
XCMA0026 
XCMA002 7 
XCMA0028 
XCMA0029 
XCMA0030 
XCMA003 1 
XCMA0032 
XCMA0033 
XCMA0034 
XCMA0035 
XCMA0036 
XCMA003 7 
XCMA0038 
XCMA0039 
XCMA0040 
XCMA004 1 
XCMA0042 
XCMA004I- 
XCMA0049 
XCMA004f 



13 CONTINUE 

BIJ=ABS(A(I »J) ) 

IF (BM-BIJ) 2 1 1 23 1 23 

21 BM = B I J 

22 IF (BM.GT.ALIMIT) GO TO 128 

23 CONTINUE . • 

I MAX= I MX 

JMAX= JMX 
AMAX=A( IMX»JMX) 

100 CONTINUE 

IF (NMAX.EO.O) GO TO 110 
•PRINT 1 

DO 121 KK=1 , NMAX 

PRINT 2t MAXI ( KK ) » MAX J ( KK ) » MAX V ( KK ) 
121 CONTINUE 
110 CONTINUE 
RETURN 
END 


XCMA0046 
XCMA004 7 
XCMA0048 
XCMA0049 
XCMA0050 
XCMA005 1 
XCMA0052 
XCMA005 3 
XCMA0054 
XCMA005 5 
XCMA0056 
XCMA005 7 
XCMA0058 
XCMA005 9 
XCMA0060 
XCMA006 1 
XCMA0062 
XCMA006 3 
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c 

c 

c 

c 

c 

c 

c 


SUBROUTINE X I JMAX ( MM , NN , I MAX , JMAX , I BM , JBM, AMAX , MAXA ) 

INPUTS TO THIS ROUTINE ARE COORDINATES AND MAGNITUDE OF PEAK VALUE 
AS FOUND IN THE WINDOW AREA CROSS-CORRE LAT I UN MATRIX XC VP ( MPT , NPT ) 
WHEN IT WAS SEARCHED BY SUBROUTINE XCMAX 

OUTPUTS OF ROUTINE ARE LAG COORDINATES QF PEAK VALUE AND 
MAGNITUDE OF CROSS-CORRELATION EXPRESSED IN PERCENTAGE UNITS 


XI JMOOO 
XI JMOOO 
XI JMOOO 
XI JMOOO 
XI JMOOO 
XI JMOOO. 
XI JMOOO 
XI JMOOO 

1 FORMAT ( 54H0 ONE PAIR VALUE IN CROSS-CORRELATION APPEARS AT I L AG= , X I JMOOO 

114, 8H JLAG= , 14 , 18H WITH PERCENT AG E=»I5»37H» WHICH WAS SELECTEDX I JMOO K 

2 AS THE BEST FIT. ) XIJMOOl 

2 FORMAT (106H0**** WARNING PEAK CROSS CORRELATION VALUE WAS FOUND AXIJMOOl. 

IT EDGE OF CORRELATION ARRAY SO MAY BE INVALID ***** ) XIJMOOl 

I BM= I MAX — 1 - ( MM/2 ) XIJMOOl' 

JBM= JMAX — 1- ( NN/2 ) XIJMOOl 

IMA=IABS( IBM ) . XIJMOOl 

JMA= I ABS ( JBM ) XIJMOOl 

IF ({ IMA. GE. (MM/2) ). OR. (JMA.GE. (NN/2) )) PRINT 2 XIJMOOl; 

MAXA=AMAX*400.0 ’ XIJMOOl' 

PRINT 1 , IBM » JBM »MAXA ■ . XIJM002< 

RETURN XIJM002 

XIJM002. 



\ 
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c 

c 

c 

c 

c 


SUBROUTINE ZFRAME ( KARA ,MM , NN , LARA , M, N, MI NCLD , MAXCLD , NF RQA , PCNTA , ZFRAOOOl 
1AMEN,SDEVA,VARA, INIT IX,TARA,NDEV) ZFRAOOOL 

DIMENSION KARA(MM,NN) ,LARA(M,N) , V AR A ( M , N ) , T AR A ( MM , NN ) ZFRA0003 

ACCEPTS AN MM BY NN ARRAY FOR TIME T-ZERO AND IDENTIFIES CENTER OUARZFR AOOO^ 
TER. THEN IT PROCESSES THAT SUBARRAY. IT NORMALIZES THE WHOLE ARKA YZFRA0005 
BY SUBTRACTING MEAN FROM CENTRAL PART AND" REPLACING FRAMING PART BY ZFRAOOUt 
ZEROES. NOTES STANDARD DEVIATION OF THE CENTRAL PART AND USES THAT ZFRA0007 
VALUE AS STANDARD DEVIATION FOR TIME T-Z ERO. M=MM/2 , N=NN/2ZFRAOOOfc 


MH=MM/4 
NH=NN/4 
DO 10 1=1, M 
DO 10 J= 1 , N 
I L = MH+ I 
JL=NH+ J 

LARA( I , J ) =KARA ( ILtJL) 

10 CONTINUE * 

CALL THRESH < L AR A , M , N ,M I NCLD , MAXClD ,NFRQA , PCNTA , AMEN , SDEVA , VARA , 
1INITIX, NDEV ) 

DO 20 1=1, MM 
DO 20 J= 1 , NN 
\TARA(I,J)=0 
20 CONTINUE 

D0\ 30 I =1 ,M \ - " 

DO 30 J=1,N \ 

I T = I +MH ^ 

JT = J + NH V 

TARA( IT ,JT)=VARA( I,J) -AMEN 
30 CONTINUE 

RETURN ' 

END 


Z FRAOOOS 
ZFRA001C 
ZFRA0011 
ZFRAOO 1 'c 
ZFRA001I- 
ZFRAOO 1 1 
ZFRA0015 
ZFRAOO 1 ( 
ZFRAOO 17 
ZFRAOOH 
ZFRAOOH 
ZFRA002C 
ZFRA002 3 
ZFRA002; 
ZFRA002: 
ZFRA002' 
ZFRA0025 
ZFRA002 ( 
ZFRA0027 
ZFRA002! 
ZFRA002‘ 
ZFRA003C 
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c 

c 

c 


SUBROUTINE ZNOKM ( KARB , MM , NN , M I NCL D , MAXCLD , NFROB , PCNTB , BMEN, 
1INITIX,TARB,NDEV) 

DIMENSION KARB(MM,NN),TARB(MM,NN) 

ACCEPTS MM BY NN ARRAY FOR TIME T-ONE , PROCESSES WITH THRESHOLDS 
M1NCLD AND MAXCLD, NORMALIZES ENTIRE ARRAY BY SUBTRACTING MEAN FROM 
EACH VALUE, NOTES STANDARD DEVIATION AND KEEPS. IT FOR TIME T/ONE 


SDEVB, ZNOROOO 1 
ZN0R0002 
ZNOR0003 
ZNOR0004 
ZN0R0005 
ZNOR0006 


CALL THRESH { K ARB , MM , NN , M I NC L D , M AX CL D, N FROB , PCNT B , BMEN , SDE VB , T ARB, ZNOROOO 7 


10 


1 1 N 1 T I X , NOE V ) 

DO 10 1=1, MM 
DO 10 J = 1 , NN 

TARBI I, J ) =TARB ( I, J) -BMEN 

CONTINUE 

RETURN 

END 


ZN0R0008 
ZN0R0009 
ZNOROO 10 
ZN0R001 1 
ZNOROO 12 
ZNOROO 1 3 
ZNOROO 14 
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Appendix G 
TESTS PERFORMED 

As described in Section 5, tests were made for the 32 array pairs 
and five control array pairs with lower threshold value of 190 and 
upper threshold values of 263, 268, and 273. Typical pages of output 
for such runs where array output printing was to be bypassed are shown 
in Figures G-l through G-5. Appendix H presents the display produced 
when the array output printing is not bypassed. 

A comprehensive timing run, whose results are summarized in 
Section 5, was conducted using Case III data. Typical pages of output 
from that run appear as Figures G-6 through G-8. 

N. • — - . * 

\ \ 



G-l 
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TAPE REAOING AND ARRAY FORMATION TOOK 0 . 05 0 0 0 SECONDS 











LATIT(JOE = 50.00000N LCNGITUOE= 6S.OOOOOW E-W MESH SIZE= 0.250000EG N- S MESH SIZE= 0.136250EG TIME BETWEEN FRAMES., 10H.MINUTFS 






J43 CMXC MAtN PROGRAM test USING PASSES 417 ANO A 1 8 NIMBUS IV T-tIR DATA 




FOLLOWING ARE TOTALS FCR ALL ARRAY PAIRS IN THIS 






Appendix H 


SAMPLE OUTPUT RESULTS 

Figures H-l through H-7 show the types of output which are obtained 
when the printing of array output is not bypassed. 
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DATA US £0 AS WINDOW AO PAV WITH COORDINATE AXES MOVED SO COORDINATES AP03AP AS LACS 



Figure H-3. WINDOW AREA DISPLAY FOR CASE VI, ARRAY PAIR 5, UPPER THRESHOLD 



DATA CONTAINED IN THE SEARCH AREA SUBARRAY IDENTIFIED AS THE BEST MATCH. COORDINATE AXFS NOVFO SO COORDINATES APPEAR AS LAGS 



Figure H-4. SEARCH AREA SUBARRAY DISPLAY FOR CASE VI. ARRAY PAIR 5. UPPER THRESHOLD - 273 



KEANS OF SUOARRAVS OF SEARCH AREA CORRESPONDING TO THF LAG VALUES 0* WINDOW AOf A PFLATfVF TO SEARCH ARFA 



Figure H-5. SUBARRAY MEAN DISPLAY FOR CASE VI, ARRAY PAIR 5, UPPER THRESHOLD 



STANDARD DEVIATION ITENTHS) OF SEARCH AREA SUBARRAY CORRESPONDIN'! TO LA?! VALUES OP WINDOW AREA RELATIVE Tn SEARCH AR p A 



Figure H-6. SUBARRAY STANDARD DEVIATION DISPLAY FOR CASE VI, ARRAY PAIR 5, UPPER THRESHOLD 





Appendix I 


SOME SUPPORTING THEORY 

1. The Cross -Co variance of Two Discrete Functions (from Stall ard) 

For two discrete functions x and y of one variable t, the cross- 
covariance function at lag x is given by 

N-l 

(1) C (t) = 2 X (t) y (t+T) 

t=o • 

where N is the number of discrete values of X and Y. This definition of 
cross-covariance is different from the modified cross -covariance defined 

cm - E 

t+0 , 

in that the unmodified form assumes that the data is cyclic, whereas the 
modified form does not. In the modified form the sum is defined only 
over those values of X and Y that overlap. This difference should be kept 
in mind when comparing the results obtained from the two different 
expressions. 

Equation (1) can be written as 

N-l N-l 

c CO = ^ ^2 s U) Y (s) 6 n ( s - t_T ) 

t=0 s=0 

where <5^ is the Kronecker delta function with its argument being considered 
modulo N, i.e. 

{ 1, if k is an integer 
0, otherwise 
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The orthogonality condition states that 

N-1 2m tu 
N 


so that 


= N 6 n (u) 


t=0 


T "4 fcl ^ 2,ir(s-t-x) 

c < t > ■ 3 .X X. x ‘ t > Y «*> X e N 

" t=0 s=0 r=0 


N-1 


N-1 


N-1 


■SUE «*>• 


-2-rri rt 
N 


»E'* < s)e 


-2m* rs \ „-2-rrirT 
e it 


r=0 


t=0 


s=0 


where Y* (S) denotes the complex conjugate of Y(S). 

Since the expressions above in parenthesis are the inverse Fourier 
transforms of X and Y respectively, 

N-1 

-2m* rx 


C(x ) = ^ A(r) B*(r) e 

r=0 


N 


where A(r) and B(r) are the Fourier coefficients of X(t) and Y(t) respectively. 

Then, by taking complex conjugates, 

N-1 

C(t) = A*(r) B(r) e . D (t). 

r=0 

where D(x) is the Fourier transform of A*B. 

For two discrete functions X and Y of two variables t-j and t 2 » the 
cross-covariance function at lag x^ and x 2 is given by 

M r l N 2 -l 

(2) C(x r x 2 ) = ^ x ( t i* V Y (VTr t 2 + l2^ 


where N^ is the number of discrete values of the variable t-j and N 2 is the 
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number of discrete values of the variable t^. Then equation (2) can be 
written as 


12 


N r’ 

N 2 -l 

N,-l 

N 2 -l 

E 

E 

E 

E 

tj =0 

t 2 =o 

s -j =0 

s 2 =0 


,<s n 1 n 2 ^ S 1 " t l “ T l* s 2 ~ x i~ x 2} 


This, 

in turn, reduces to 



V 1 

n 2 -! 

Vl r i T l 

C (x-| 

• l 2> = E 

T: A(r ] ,r 2 ) B* (r r r 2 ) e " 27T1 

N, N, 


r r° 

0 

11 

CM 

V. 



N i — 1 N 2 -l 


c(t r T 2> = E E A * (r i> r 2 ) B (r r r 2 ) e 



2\ . 2ui 


r ri ^ r 2 T 2 


r f0 r 2 = 0 


A. 


which is D ( t j, t 2 ), the Fourier transform of A*B. 
2. Covariance Computation (From Stall ard) 


The classical means of computing the cross-covariance function had 
been lagged products. However, the advent of a fast Fourier transform allowed 
a new approach to computing the cross-covariance. 


The cross -covariance of two functions f and g at lags u^ and u 2 is 
approximated by 

C(u r u 2 ) = T" 1 C(T ( f ) )* T(g)] 

where T is the finite Fourier series representation of a function. For the 
case in which f and g are the same functions the function C(u^, u 2 ) is 
referred to as the auto-covariance. 
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If f represents one set of cloud data, and g represents a subsequent 
set of cloud data, then the motion required to move f with respect to g so 
as to produce a maximum covariance- is the motion of the cloud pattern. 

3. Simultaneous Fourier Analysis of Two Sets of Real Data (From Cooley) 

Here we shall describe a procedure for calculating the Fourier 
transforms of two sets of real data by applying one complex Fourier trans- 
form. Letting- S^ (j) and X 2 (j) be two sets of real data with 





X-,(j) «-*. A^n) 
X 2 ( J ) •*-*■ A 2 ( n ) , 

we have 

X(j) ~ A(n) 

where 

X(j) = X-j(j) + 

and 


hence 

and 


A(n) = A-j(n) + i A 2 (n). 

Replacing n by N-n, taking conjugates of both sides we may obtain 
A (H-n) = A-j i (N-n) - i A 2 (N-n) = A-j(n) - i A 2 (n) 

It then follows by addition that 

A(n) + A (N-n) = 2A-j(n) and by subtraction 
A(n) - A*(N-n) = 2i A 2 (n) 

A, (n) = 1 [A(n) + A* (N-n)] 

A 2 (n) = [A(n) - A* (KMi)] 

For the special cases n=0 and n=N/2 we note that 
A(0) = A, (0) + 1 A 2 (0) 

A(£) = Ai (|) + i A 2 (|) 
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So A-j(O) is the Real part of A(0) 

A 2 (0) is the Imaginary part of A(0) 

A^ (^-) is the Real part of A(^-) 

N N 

A 2 is the Imaginary part of A^) 

Since X-j(j) and X 2 (j) are postulated as two sets of real data their 
Fourier coefficients satisfy the symmetry property 

A-j(n) = A-j (N-n) and A 2 (n) = (N-n) 

and therefore only one half of each array need be computed and stored. 

Thus one requires the same amount of storage for A^(n) and A^( n ) as for 
A(n) or X(j). A suggested storage arrangement is to replace A(n) by A-j(n) 
and A(N-n) by A 2 (N-n) for n = 1, 2, . .., N/2-1. We see that we can skip 
the indices n = 0 and N/2 since we already have the results in the’ real 
and imaginary parts of A(0) and A(N/2). To summarize the procedure when 
using the subroutine HARM: 

Step 1. Let X (j) = X-j(j) + i X 2 (j) be the input to HARM 

Step 2. Call HARM, requesting that a Fourier transform be computed, 

replacing the array X(.i) by A(n). 

Step 3. For n = 1, 2, ..., (M/2 -1), let 

A-j(n) = 2 - [A(n) + A*(N-n)] 

A^N-n) = [A(n) - A*(N-n}] 

with A^(n) and A 2 (N-n) replacing A(n) and A(N-n), respectively, 
in storage. 

4. Errors in Digital Fourier Transform Convolution/Correlation 
(From Silverman) 

Errors in performing a DFT convolution or correlation may be 
considered as stemming from two sources. First is the error incurred in 
sampling continuous functions, which is easy to bound. Second is the error 



resulting from roundoff in fast Fourier transforms and in multiplications. 

: j;-:;; aliasing to he smal 1 ’■?■/ ’• nclurfino the vast portico of the region of 

The error in approximation of, a continuous Fourier Transform .by the 


it :u;s ' 


:?o !ar-'H* enough so t.ia~ 


DFT has been shown to result- from three sources. First, one assumes that 

tue ' reeuence acme*.?; wi ? ; t v nn« enounn. iho parameter c •r.t.s a iso oe 

the function is non-zero only over a finite interval. Should this not be 

:> (* ■, o c te Q to opt >r.n te two erf ten a. 1 1 may ne selected to nirunrze tne 

the case, aliasinq error occurs.. The magnitude of this error can easily 

arpro;om.K8 arrorV-r tee triancusar •snteuremen , wr.? is at the sem-e f.me 

be bounded by comparing a norm for the function on the intervals, (-» ,0), 

’nsunna W3 u cne ranoe un toe frep’.'oncv domain is jarcc. . Or course, tnera 

(T, oo) with the same norm over (0,T). This is demonstrated for the square 

f > r cost assoc* a .eu wun maxmq : and sma < ,* in terms of comouta- 


norm in 


-V / SC'HS to {'i'Vs r - i d' : * .'i '«yVr\ V". v> r t ii tot "i £ ! f* d i 




" llx(t)! ! ab d t7 X 2 (t)dt 3 ^ ' • • 

The errors due uo rou:i 7 'ff in calculation of the correlation 


or con VO 


! Ul w " % aliasinq error < 100 x ’( 1 1 X( t) ] 1 " ‘ n + 1 1 X(t) 1 1 ’ ) ; 

vf one perror-ns jj scrota corra iac-o n i1v ' t 

plicati vs process can yield an accumulati vn lr.Q.,:T;rti onal to the size 

;‘(N.o.te:u? there are sped al .cases. when rthe ; regi on /-r,» , °o ) need .not .be 
considered.) . . , • . 

y.f' £.r\p Y'OO V 1 d V.'rl i V 11 'l f'UCt'u rOO Or ; D i A 0 Cr 3 tTi ’.Vnl Ch CO 0 S OOt :*!U , !Vi 0 i V 

for A second : aliasing error occurs; as a resul t of saving only a -fini te 
n.umber.iOf frequencies, o If an tanalyjtic ; solution /or- the. continuous /purler 
■Transform ,for -a,. /unction, is .known .ythe -application of /he.. above in equality 
/or .frequency rdomain. /unctions-'.mayi/e appl ied ftp, find this : , source vof error 
■by. c.repl acing i'T/;Wi.th /vans form. This general ‘averaging pronerty of the 
truncr^. ^ r d source of'iefrdr iif r thiV l appfox'i matidrf comes' c ffom r the‘ c trape-’ 
r zoi dal 'i ntegra ti on , uS /hints' Wop^ 'the Second deri vati ve of D 1 ‘ " 

‘the '-'function being approximated. 

Experience has shown ths l truncation error does not general ly accumulate 


U\3 


712 


- /" (nAt) : . orru * 5 ci ens/ convc ; ufi ans when r * cad ng-po ! n v. 


arithmetic ‘is used/ There seems no point in a voiding the use of DFT -solution 


ThereforeY 'the DFT is r a reasonable ' approximatidh to the continuous Fourier 


Transform only if all three of these sources of error are small. These 
three error criteria severely influence the selection of the two parameters 
which may be selected, t and T. The integration interval T must be 
selected to optimize two criteria. First, it must cause the error due to 
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